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Inspirals of stellar mass compact objects into massive black holes are an important source for 
future gravitational wave detectors such as Advanced LIGO and LISA. Detection of these sources and 
extracting information from the signal relies on accurate theoretical models of the binary dynamics. 
We cast the equations describing binary inspiral in the extreme mass ratio limit in terms of action 
angle variables, and derive properties of general solutions using a two-timescale expansion. This 
provides a rigorous derivation of the prescription for computing the leading order orbital motion. As 
shown by Mino, this leading order or adiabatic motion requires only knowledge of the orbit-averaged, 
dissipative piece of the self force. The two timescale method also gives a framework for calculating 
the post-adiabatic corrections. For circular and for equatorial orbits, the leading order corrections 
are suppressed by one power of the mass ratio, and give rise to phase errors of order unity over a 
complete inspiral through the relativistic regime. These post-l-adiabatic corrections are generated 
by the fluctuating, dissipative piece of the first order self force, by the conservative piece of the first 
order self force, and by the orbit-averaged, dissipative piece of the second order self force. We also 
sketch a two-timescale expansion of the Einstein equation, and deduce an analytic formula for the 
leading order, adiabatic gravitational waveforms generated by an inspiral. 



I. INTRODUCTION AND SUMMARY 
A. Background and Motivation 

Recent years have seen great progress in our under- 
standing of the two body problem in general relativ- 
ity. Binary systems of compact bodies undergo an inspi- 
ral driven by gravitational radiation reaction until they 
merge. As illustrated in Fig. [U there are three different 
regimes in the dynamics of these systems, depending on 
the values of the total and reduced masses M and ^ of the 
system and the orbital separation r : (i) The early, weak 
field regime at r ^ M, which can be accurately modeled 
using post-Newtonian theory, see, for example, the re- 
view [l|. (ii) The relativistic, equal mass regime r ~ M, 
fi ~ M, which must be treated using numerical relativ- 
ity. Over the last few years, numerical relativists have 
succeeded for the first time in simulating the merger of 
black hole binaries, see, for example, the review 2-1 and 
references therein, (iii) The relativistic, extreme mass 
ratio regime r ^ M, fi M. Over times cales short com- 
pared to the dephasing time '-^ M^/M/fj,, systems in this 
regime can be accurately modeled using black hole per- 
turbation theory [3!], with the mass ratio e = n/M serving 
as the expansion parameter. The subject of this paper is 
the approximation methods that are necessary to treat 
such systems over the longer inspiral timescale ~ M^//i 
necessary for computation of complete inspirals. 

This extreme mass ratio regime has direct observa- 
tional relevance: Compact objects spiraling into much 
larger black holes are expected to be a key source for 
both LIGO and LISA. Intermediate-mass-ratio inspirals 
(IMRIs) are inspirals of black holes or neutron stars into 
intermediate mass (50 < M < IOOOMq) black holes; 
these would be visible to Advanced LIGO out to dis- 
tances of several hundred Mpc , where the event rate 
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FIG. 1: The parameter space of inspiralling compact binaries 
in general relativity, in terms of the inverse mass ratio M/ fi = 
1/e and the orbital radius r, showing the different regimes 
and the computational techniques necessary in each regime. 
Individual binaries evolve downwards in the diagram (green 
dashed arrows). 



could be about 3 — 30 per year d, [^. Extreme-mass- 
ratio inspirals (EMRIs) are inspirals of stellar-mass com- 
pact objects (black holes, neutron stars, or possibly white 
dwarfs) into massive (lO'' < M < lO^M©) black holes 
in galactic nuclei; these will be visible to LISA out to 
redshifts z w 1 [1, 0, 1!]. It has been estimated [1 [l3| 
that LISA should see about 50 such events per year, 
based on calculations of stellar dynamics in galaxies' cen- 
tral cusps [llj. Because of an IMRFs or EMRFs small 
mass ratio e — n/M, the small body lingers in the large 
black hole's strong-curvature region for many wave cy- 
cles before merger: hundreds of cycles for LIGO's IM- 
RIs; hundreds of thousands for EISA's EMRIs 0. In 
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this relativistic regime the post-Newtonian approxima- 
tion has completely broken down, and full numerical rel- 
ativity simulations become prohibitively difficult as e is 
decreased. Modeling of these sources therefore requires 
a specialized approximation method. 

Gravitational waves from these sources will be rich 
with information 0, : 

• The waves carry not only the details of the evolving 
orbit, but also a map of the large body's spacetime 
geometry, or equivalently the values of all its mul- 
tipole moments, as well as details of the response 
of the horizon to tidal forces [H, ■ Extracting 
the map (bothrodesy) is a high priority for LISA, 
which can achieve ultrahigh accuracy, and a mod- 
erate priority for LIGO, which will have a lower 
(but still interesting) accuracy y. Measurements 
of the black hole's quadrupole (fractional accuracy 
about 10-3 for LISA [li,[ll], about 1 for Advanced 
LIGO ^^l) will enable tests of the black hole's no 
hair property, that all of the mass and current mul- 
tipole moments are uniquely determined in terms 
of the first two, the mass and spin. Potentially, 
these measurements could lead to the discovery of 
non-black-hole central objects such as boson stars 

, [l3| or naked singularities. 

• One can measure the mass and spin of the cen- 
tral black hole with fractional accuracies of order 
10--* for LISA and about IQ-^-lQ-^ for 
Advanced LIGO [41]. Observing many events will 
therefore provide a census of the masses and spins 
of the massive central black holes in non-active 
galactic nuclei like M31 and M32. The spin can 
provide useful information about the hole's growth 
history (mergers versus accretion) (20l |. 

• For LISA, one can measure the inspiralling ob- 
jects' masses with precision about 10^^, teaching 
us about the stellar population in the central par- 
sec of galactic nuclei. 

• If the LISA event rate is large enough, one can 



measure the Hubble constant Hq to about 1% 21 1, 
which would indirectly aid dark energy studies [22 1 . 
The idea is to combine the measured luminosity 
distance of cosmological {z ^ 1/2) EMRIs with a 
statistical analysis of the redshifts of candidate host 
galaxies located within the error box on the sky. 



To realize the science goals for these sources requires 
accurate theoretical models of the waveforms for matched 
filtering. The accuracy requirement is roughly that the 
theoretical template's phase must remain accurate to ^ 1 
cycle over the cycles of waveform in the highly 

relativistic regime (~ 10^ cycles for LIGO, ^ 10^ for 
LISA). For signal detection, the requirement is slightly 
less stringent than this, while for parameter extraction 
the requirement is slightly more stringent: The wave- 
forms are characterized by 14 parameters, which makes 



a fully coherent search of the entire data train compu- 
tationally impossible. Therefore, detection templates for 
LISA will use short segments of the signal and require 
phase coherence for ~ 10'' cycles [T^]. Once the presence 
of a signal has been established, the source parameters 
will be extracted using measurement templates that re- 
quire a fractional phase accuracy of order the reciprocal 
of the signal to noise ratio [l^l , in order to keep system- 
atic errors as small as the statistical errors. 



B. Methods of computing orbital motion and 
waveforms 



A variety of approaches to computing waveforms have 
been pursued in the community. We now review these ap- 
proaches in order to place the present paper in context. 
The foundation for all approaches is the fact that, since 
e = iijM <^ 1, the field of the compact object can be 
treated as a small perturbation to the large black hole's 
gravitational field. On short timescales ~ M, the com- 
pact object moves on a geodesic of the Kerr geometry, 
characterized by its conserved energy E, z-component of 
angular momentum L^, and Carter constant Q. Over 
longer timescales ~ M/e, radiation reaction causes the 
parameters E, and Q to evolve adiabatically and the 
orbit to shrink. The effect of the internal structure of 
the object is negligible^, so it can be treated as a point 
particle. At the end of the inspiral, the particle passes 
through an innermost stable orbit where adiabaticity 
breaks down, and it transitions onto a geodesic plunge 
orbit [13, m, [2^ [13] ■ In this paper we restrict attention 
to the adiabatic portion of the motion. 

Numerical Relativity: Numerical relativity has not yet 
been applied to the extreme mass ratio regime. How- 
ever, given the recent successful simulations in the equal 
mass regime e ~ 1, one could contemplate trying to per- 
form simulations with smaller mass ratios. There are a 
number of difficulties that arise as e gets small: (i) The 
orbital timescale and the radiation reaction timescale are 
separated by the large factor ^ 1/e. The huge number of 
wave cycles implies an impractically large computation 
time, (ii) There is a separation of lengthscales: the com- 
pact object is smaller than the central black hole by a 
factor £. (iii) Most importantly, in the strong field region 
near the small object, the piece of the metric perturba- 
tion responsible for radiation reaction is of order e, and 



^ There are two exceptions, where corrections to the point-particle 
model can be important: (i) White dwarf EMRIs, where tidal 
interactions can play a role [23II . (ii) The effect due to the 
spin, if any, of the inspirallin g ob ject, whose importance has 
been emphasized by Burko l24l |25|| . While this effect is at most 
marginally relevant for signal detection |26l| , it is likely quite im- 
portant for information extraction. We neglect the spin effect in 
the present paper, since it can be computed and included in the 
waveforms relatively easily. 
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since one requires errors in the radiation reaction to be of 
order e, the accuracy requirement on the metric pertur- 
bation is of order . These difficulties imply that numer- 
ical simulations will likely not be possible in the extreme 
mass ratio regime in the foreseeable future, unless major 
new techniques are devised to speed up computations. 

Use of post-Newtonian methods: Approximate waveforms 
which are qualitatively similar to real waveforms can 
be obtained using post-Newtonian methods or using hy- 
brid schemes containing some post-Newtonian elements 
[H, El, mi ■ Although these waveforms are insufficiently 
accurate for the eventual detection and data analysis of 
real signals, they have been very useful for approximately 
scoping out the detectability of inspiral events and the 
accuracy of parameter measurement, both for LIGO Q 
and LISA [13, [2^ . They have the advantage that they 
can be computed relatively quickly. 

Black hole perturbation theory - first order: There is 
a long history of using first order perturbation theory 
to compute gravitational waveforms from par ticles 
in geodesic orbits around black holes [H, [13, ISSl . [36j . 
These computations have recently been extended to fully 
generic orbits [13, [H, [s^. However first order pertur- 
bation theory is limited to producing "snapshot" wave- 
forms that neglect radiation reaction.^ Such waveforms 
fall out of phase with true waveforms after a dephasing 
time ~ M/y/e, the geometric mean of the orbital and ra- 
diation reaction timescales, and so are of limited utility.^ 

Black hole perturbation theory - second order: One can 
in p rinciple go to second order in perturbation theory 
|4ll. |4^. [43} . At this order, the particle's geodesic mo- 
tion must be corrected by self-force effects describing 
its interaction with its own spacetime distortion. This 
gravitational self force is analogous to the electromag- 
netic Abraham-Lorentz-Dirac force. Although a formal 
expression for the self force is known [3, |45| | , it has 
proved difficult to translate this expression into a practi- 
cal computational scheme for Kerr black holes because 
of the mathematical complexity of the self-field regu- 
larization which is required. Research into this topic 
is ongo ing ; see, for example the review [i^ and Refs. 
m, [43, HE m, [m, [is, ^r various approaches 

and recent progress. 

When the self-force is finally successfully computed, 
second order perturbation theory will provide a self- 
consistent framework for computing the orbital motion 



and the waveform, but only over short timescales. The 
second order waveforms will fall out of phase with the 
true waveforms after only a dephasing time ^ Mj^ ^ 
[55 ^, '5^ . Computing accurate waveforms describing a full 
inspiral therefore requires going beyond black hole per- 
turbation theory. 

Use of conservation laws: This well-explored method al- 
lows tracking an entire inspiral for certain special classes 
of orbits. Perturbation theory is used to compute the 
fluxes of E and to infinity and down the horizon for 
geodesic orbits, and imposing global conservation laws, 
one infers the rates of change of the orbital energy and 
angular momentum. For circular orbits and equatorial 
orbits these determine the rate of change of the Carter 
constant Q, and thus the inspiralling trajectory. The 
computation can either be done in the frequency domain 
[33I [33, [3^ , or in the time domain by numerically 
integrating the Teukolsky equation as a 2+1 PDE with 
a suitable numerical model of the point particle source 
[13, [H, Hi, [13, HH M, M, Hi, [ei, H^. However, this 
method fails for generic orbits since there is no known 
global conservation law associated with the Carter con- 
stant Q. 

Adiabatic approximation - leading order: Over the last 
few years, it has been discovered how to compute inspi- 
rals to leading order for generic orbits. The method is 
based on the Mino's realization [g^] that, in the adia- 
batic limit, one needs only the time averaged, dissipative 
piece of the first order self force, which can be straight- 
forwardly computed from the half retarded minus half 
advanced prescription. This sidesteps the difficulties as- 
sociated with regularization that impede computations 
of the full, first order self force. From the half advanced 
minus half retarded prescription, one can derive an ex- 
plicit formula for a time-average of Q in terms of mode 
expansion [13, [H, [13, [tO, [71| . Using this formula it will 
be straightforward to compute inspirals to the leading 
order. 

We now recap and assess the status of these various 
approaches. All of the approaches described above have 
shortcomings and limitations [56| . Suppose that one 
computes the inspiral motion, either from conservation 
laws, or from the time-averaged dissipative piece of the 
first order self-force, or from the exact first order self- 
force when it becomes available. It is then necessary to 
compute the radiation generated by this inspiral. One 
might be tempted to use linearized perturbation theory 
for this purpose. However, two problems then arise: 



^ The source for the hnearized Einstein equation must be a con- 
served stress energy tensor, which for a point particle requires a 
geodesic orbit. 

^ Drasco has argued that snapshot waveforms may still be useful 
for signal detections in certain limited parts of the IMRI/EMRI 
parameter space, since the phase coherence time is actually ^ 
lOOM/v^ [3. 



^ The reason is as follows. Geodesic orbits and true orbits become 
out of phase by ~ 1 cycle after a dephasing time. Therefore, since 
the linear metric perturbation is sourced by a geodesic orbit, 
fractional errors in the linear metric perturbation must be of 
order unity. Therefore the second order metric perturbation must 
become comparable to the first order term after a dephasing time. 
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• As noted above, the use of linearized perturba- 
tion theory with nongeodesic sources is mathemat- 
ically inconsistent. This inconsistency has often 
been remarked upon, and various ad hoc methods 
of modifying the linearized theory to get around 
the difficulty have been suggested or implemented 

[ii,[z3,[zi. 

• A related problem is that the resulting waveforms 
will depend on the gauge chosen for the linearized 
metric perturbation, whereas the exact waveforms 
must be gauge invariant. 

It has often been suggested that these problems can be 
resolved by going to second order in perturbation theory 
(43I l46j . However, as discussed above, a second order 
computation will be valid only for a dephasing time, and 
not for a full inspiral. 

Of course, the above problems are not fatal, since the 
motion is locally very nearly geodesic, and so the incon- 
sistencies and ambiguities are suppressed by a factor ~ e 
relative to the leading order waveforms.^ Nevertheless, 
it is clearly desirable to have a well defined approxima- 
tion method that gives a unique, consistent result for the 
leading order waveform. Also, for parameter extraction, 
it will be necessary to compute the phase of the waveform 
beyond the leading order. For this purpose it will clearly 
be necessary to have a more fundamental computational 
framework. 



C. The two timescale expansion method 

In this paper we describe an approximation scheme 
which addresses and resolves all of the theoretical dif- 
ficulties discussed above. It is based on the fact that 
the systems evolve adiabatically: the radiation reac- 
tion timescale ^ M/e is much longer than the orbital 
timescale ~ M [g^]- It uses two-timescale expansions, 
which are a systematic method for studying the cumula- 
tive effect of a small disturbance on a dynamical system 
that is active over a long time (t^ . 

The essence of the method is simply an ansatz for the 
dependence of the metric gabis) on e, and an ansatz for 
the dependence of the orbital motion on e, that are jus- 
tified a posteriori order by order via substitution into 
Einstein's equation. The ansatz for the metric is more 
complex than the Taylor series ansatz which underlies 
standard perturbation theory. The two timescale method 
has roughly the same relationship to black hole pertur- 
bation theory as post-Newtonian theory has to pertur- 
bation theory of Minkowski spacetime. The method is 
consistent with standard black hole perturbation theory 
locally in time, at each instant, but extends the domain 



^ This is true both for the instantaneous amphtude and for the 
accumulated phase of the waveform. 



of validity to an entire inspiral. The method provides 
a systematic procedure for computing the leading order 
waveforms, which we call the adiabatic waveforms, as 
well as higher order corrections. We call the 0(e) cor- 
rections the post-l-adiabatic corrections, the 0(£^) cor- 
rections post-2-adiabatic, etc., paralleling the standard 
terminology in post-Newtonian theory. 

The use of two timescale expansions in the extreme 
mass ratio regime was first suggested in Refs. [H, |75| . 
The method has already been applied to some simplified 
model problems: a computation of the inspiral of a point 
particle in Schwarzschild subject to electromagnetic ra- 
diation reaction forces by Pound and Poisson [t^, and 
a computation of the scalar radiation generated by a in- 
spiralling particle coupled to a scalar field by Mino and 
Price [73]. We will extend and generalize these analyses, 
and develop a complete approximation scheme. 

There are two, independent, parts to the the approx- 
imation scheme. The first is a two timescale analysis 
of the inspiralling orbital motion, which is the focus of 
the present paper. Our formalism enables us to give a 
rigorous derivation and clarification of the prescription 
for computing the leading order motion that is valid for 
all orbits, and resolves some controversies in the litera- 
ture [7^. It also allows us to systematically calculate 
the higher order corrections. For these corrections, we 
restrict attention to inspirals in Schwarzschild, and to 
circular and equatorial inspirals in Kerr. Fully generic 
inspirals in Kerr involve a qualitatively new feature - the 
occurrence of transient resonances - which we will discuss 
in the forthcoming papers [78, 79] . 

The second part to the approximation scheme is the 
application of the two timescale method to the Einstein 
equation, and a meshing of that expansion to the analysis 
of the orbital motion. This allows computation of the 
observable gravitational waveforms, and is described in 
detail in the forthcoming paper [8^. We briefly sketch 
this formalism in Sec. II El below, and give an analytic 
result for the leading order waveforms. 

We note that alternative methods of attempting to 
overcome the problems with standard perturbation the- 
ory, and of going beyond adiabatic order, have been de- 
veloped by Mino [11 [zl, [8l|, HI, HI- These methods 
have some overlap with the method discussed here, but 
differ in some crucial aspects. We do not believe that 
these methods provide a systematic framework for going 
to higher orders, unlike the two-timescale method. 



D. Orbital Motion 

We now turn to a description of our two timescale anal- 
ysis of the orbital motion. The first step is to exploit the 
Hamiltonian structure of the unperturbed, geodesic mo- 
tion to rewrite the governing equations in terms of gen- 
eralized action angle variables. We start from the forced 
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(i)- + £2a(2)- + 0(e3). (1.1) 



geodesic equation 
(fx" dx" dxP _ 

Here r is proper time and a*^^-' and a*^^^ " are the first 
order and second order self-accelerations. In Sec. [Til we 
augment these equations to describe the leading order 
backreaction of the inspiral on the mass M and spin a 
of the black hole, and show they can be rewritten as [cf. 
Eqs. (|2Tf|l below] 

+0(£^), (1.2a) 

eG^^\qr,q0, J a) + £^G^^'(gr, qe, J a) 
+0(£3). (1.2b) 



dJx 
dT 



Here the variables J\ are the three conserved quantities 
of geodesic motion, with the dependence on the particle 
mass scaled out, together with the black hole mass and 
spin parameters: 



Jx = {E/fi,L,/fi,Q/fi^,M,a). 



(1.3) 



The variables qa — {qr,q9,q4>,qt) are a set of general- 
ized angle variables associated with the r, 9, (j) and t 
motions in Boyer-Lindquist coordinates, and are defined 
more precisely in Sec. IIIDI below. The variables qr, qe, 
and q^ each increase by 27r after one cycle of motion of the 
corresponding variable r, 6 or (jj. The functions LUa{Jcr) 
are the fundamental frequencies of geodesic motion in the 
Kerr metric. The functions ga \ G^^^ are currently not 
known explicitly, but their exact form will not be needed 
for the analysis of this paper. They arc determined by 
the first order self acceleration [4j, |45|. Similarly, the 

(2) (2) 

functions ga and G\ are currently not known explic- 
itly, and are determined in part by the second order self 
acceleration [H HE [H, [83, ill; see Sec. [HFl for more 
details. 

In Sees. llVI -rVlbelow we analyze the differential equa- 
tions (|1.2p using two timescale expansions. In the non- 
resonant case, and up to post-l-adiabatic order, the re- 
sults can be summarized as follows. Approximate solu- 
tions of the equations can be constructed via a series of 
steps: 

• We define the slow time variable f = er. 

• We construct a set of functions Tpa\f), Jx^\f), 

ipa\f) and Jx^\f) of the slow time. These func- 
tions are defined by a set of differential equations 



ga^ and 



that involve the functions uja, ga \ G^^\ 

G\ and which are independent of e [Eqs. (|5.26p . 
Km . Km . (|539)1 . (lOTll below]. 

• We define a set of auxiliary phase variables ipa by 



where the 0{e) symbol refers to the limit e 
fixed f = ST. 



at 



• Finally, the solution to post-l-adiabatic order is 
given by 



^pa+Oie), (1.5a) 
+Hx[A, Ve, J-j"^ (er)] + 0{e^), (1.5b) 



where the 0{e) and O(e^) symbols refer to e ^ 
at fixed f and ipa- Here H\ is a function which is 
periodic in its first two arguments and which can 
computed from the function G^^'' [Eq. (17. 3|) below]. 

We now turn to a discussion of the implications of 
the final result (jl.5p . First, we emphasize that the pur- 
pose of the analysis is not to give a convenient, practical 
scheme to integrate the orbital equations of motion. Such 
a scheme is not needed, since once the self-acceleration is 
known, it is straightforward to numerically integrate the 
forced geodesic equations (|l.ip . Rather, the main benefit 
of the analysis is to give an analytic understanding of the 
dependence of the motion on e in the limit e — > 0. This 
serves two purposes. First, it acts as a foundation for the 
two timescale expansion of the Einstein equation and the 
computation of waveforms (Sec. II El below and Ref. [80|). 
Second, it clarifies the utility of different approximations 
to the self-force that have been proposed, by determining 
which pieces of the self-force contribute to the adiabatic 
order and post-l-adiabatic order motions [s^. 111]. This 
second issue is discussed in detail in Sec. lVIII below. Here 
we give a brief summary. 

Consider first the motion to adiabatic order, given by 
the functions ipa ^ and J'^^^ . These functions are obtained 
from the differential equations [Eqs. (|5.26p . (|5.3ip and 
([Og)) below] 



dT 



(1.6a) 



dj. 



(0) 



df 



ir) = {Gi'^)[J^'\f)], (1.6b) 



where 



denotes the average^ over the 2-torus 



(2^)2 



dqr 



dqe G^^\qr,qe,Ja)- 



(1-7) 

This zeroth order approximation describes the inspi- 
ralling motion of the particle. In Sec. Ill Gl below we show 
that only the dissipative (ie half retarded minus half ad- 
vanced) piece of the self force contributes to the torus 



® This phase space average is uniquely determined by the dynam- 
ics of the system, and resolves concerns in the literature about 
inherent ambiguities in the choice of averaging |76| . 
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average (11.71) . Thus, the leading order motion depends 
only on the dissipative self- force, as argued by Mino [67j . 
Our result extends slightly that of Mino, since he advo- 
cated using an infinite time average on the right hand 
side of Eq. (|1.6bp . instead of the phase space or torus 
average. The two averaging methods are equivalent for 
generic geodesies, but not for geodesies for which the 
ratio of radial and azimuthal periods is a rational num- 
ber. The time-average prescription is therefore correct 
for generic geodesies, while the result (|1.6p is valid for 
all geodesies. The effect of the nongeneric geodesies is 
discussed in detail in Refs. [tsLItqI. 

Consider next the subleading, post-l-adiabatic correc- 
tions to the inspiral given by the functions ipa^ and J^^^ ■ 
These corrections are important for assessing the accu- 
racy of the adiabatic approximation, and will be needed 
for accurate data analysis of detected signals. The dif- 
ferential equations determining tpa'^ and JT^^"* are Eqs. 
(|5.39p and (|5.37p below. These equations depend on (i) 
the oscillating (not averaged) piece of the dissipative, first 
order self force; (ii) the conservative piece of the first or- 
der self force, and (iii) the torus averaged, dissipative 
piece of the second order self force. Thus, all three of 
these quantities will be required to compute the inspi- 
ral to subleadin g or der, confirming arguments made in 
Refs. jsrllssllsgnQOt . In particular, knowledge of the full 
first order self force will not enable computation of more 
accurate inspirals until the averaged, dissipative piece of 
the second order self force is known. ^ 



E. Two timescale expansion of the Einstein 
equations and adiabatic waveforms 

We now turn to a brief description of the two timescale 
expansion of the Einstein equations; more details will be 
given in the forthcoming paper js^l- We focus attention 
on a region TZ of spacetime defined by the conditions (i) 
The distance from the particle is large compared to its 
mass fi; (ii) The distance r from the large black hole is 
small compared to the inspiral time, r <C M'^/fi; and 
(iii) The extent of the region in time covers the entire 
inspiral in the relativistic regime. In this domain we make 
an ansatz for the form of the metric that is justified a 
posteriori order by order. 

At distances ~ fi from the particle, one needs to use 
a different type of analysis (eg black hole perturbation 
theory for a small black hole), and to mesh that analy- 
sis with the solution in the region TZ by matching in a 
domain of common validity. This procedure is very well 
understood and is the standard method for deriving the 
first order self force [3, [9l[ . It is valid for our metric 
ansatz (II. 8[) below since that ansatz reduces, locally in 



time at each instant, to standard black hole perturbation 
theory. Therefore we do not focus on this aspect of the 
problem here. Similarly, at large distances, one needs 
to match the solution within TZ onto an outgoing wave 
solution in order to read off the asymptotic waveforms.^ 
Within the region TZ, our ansatz for the form of the 
metric in the non-resonant case is 



+e^gl^i{qr,qe,q^,i,x3) + 0{e^). (1.8) 

Here g^^^ is the background, Kerr metric. The coordi- 
nates (t, x^ ) can be any set of coordinates in Kerr, subject 
only to the restriction that d/dt is the timelike Killing 
vector. On the right hand side, t is the slow time variable 
t = et, and the quantities qr, qe and q^ are the values of 
the orbit's angle variables at the intersection of the inspi- 
ralling orbit with the hypersurface t = constant. These 
are functions of t and of e, and can be obtained from 
the solutions p.4p and (|1.5ap of the inspiral problem by 
eliminating the proper time r. The result is of the form 



-fl°\i) + f^'\i) + 0{e), (1.9) 



for some functions f^^\ fl^\ On the right hand side of 
Eq. (|1.8p . the 0{£^) refers to an asymptotic expansion 
associated with the limit e ^ at fixed g^, x^ and i. 



Finally the functions g^^ 



and g^^p are assumed to be 



multiply periodic in g,., qe and with period 27r in each 
variable. 

By inserting the ansatz (|1.8p into Einstein's equations, 
one obtains a set of equations that determines the free 
functions, order by order. At leading order we obtain an 
equation of the form 



0, 



(1.10) 



where I? is a linear differential operator on the six di- 
mensional manifold with coordinates {qr,qe,q4>,x^)- In 
solving this equation, t is treated as a constant. The 
solution that matches appropriately onto the worldline 
source can be written as 



6M{t) 



6a{i) 



dM ' ' da 
+^ai3[qr, qe,q<j,, x\E{t), L^{t),Q(i)]. 



(1.11) 



Here the terms on the first line are the secular pieces of 
the solution. They arise since the variable t is treated as 
a constant, and so one can obtain a solution by taking 
the perturbation to the metric generated by allowing the 



This statement remains true when one takes into account reso- 
nances [79II . 



This matching is not necessary at the leading, adiabatic order, 
for certain special choices of time coordinate in the background 
spacetime, as argued in Ref. [TtII . It is needed to higher orders. 
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parameters of the black hole (mass, spin, velocity, center 
of mass location) to vary as arbitrary functions of t. For 
example, the mass of the black hole can be written as 
M{i) = M + 6M{i), where M = M{0) is the initial mass. 
The functions SM{i), Sa{i) etc. are freely specifiable at 
this order, and will be determined at the next (post-1- 
adiabatic) order. 

The second line of Eq. (jl.lip is the oscillatory piece of 
the solution. Here one obtains a solution by taking the 
function JF^^ to be the function 

^apiqr, qe,q^, x\E, L^, Q) 

that one obtains from standard linear perturbation the- 
ory with a geodesic source. This function is known 
analytically in Boyer-Lindquist coordinates (t, r, 9, 0) in 
terms of a mode expansion.^' 

The gauge freedom in this formalism consists of those 
one parameter families of diffeomorphisms which preserve 
the form (|1.8p of the metric ansatz. To the leading order, 
these consist of (i) gauge transformations of the back- 
ground coordinates that are independent of e, which pre- 
serve the timelike Killing vector, and (ii) transformations 
of the form 

^ + (^Zr, ge, 90, i, x^) + Oie^). (1.12) 

Note that this is not the standard gauge freedom of lin- 
ear perturbation theory, since depends on 4 "time 
variables" instead of one. This modified gauge group al- 
lows the two timescale method to evade the two problems 
discussed at the end of Sec. II Bl above, since the gradual 
evolution is described entirely by the t dependence, and, 
at each fixed i, the leading order dependence on the vari- 
ables qr, qg, q,p, r, 6 and ip is the same as in standard 
perturbation theory with the same gauge transformation 
properties. 



F. Organization of this Paper 

The organization of this paper is as follows. In Sec. |TT] 
we derive the fundamental equations describing the in- 
spiral of a point particle into a Kerr black hole in terms of 
generalized action-angle variables. In Sec. IIIII we define 
a class of general, weakly perturbed dynamical systems 
of which the inspiral motion in Kerr is a special case. 
We then study the solutions of this class of systems us- 
ing two-timescale expansions, first for a single degree of 
freedom in Sec. IIVI and then for the general case in Sec. 
FVl Section IVII gives an example of a numerical integra- 
tion of a system of this kind, and Sec. IVIII gives the final 
discussion and conclusions. 



G. Notation and Conventions 

Throughout this paper we use units with G = c = 1. 
Lower case Roman indices a. b, c, . . . denote abstract in- 
dices in the sense of Wald [94|. We use these indices 
both for tensors on spacetime and for tensors on the 
eight dimensional phase space. Lower case Greek indices 
ly, A, (T, T, . . . from the middle of the alphabet denote com- 
ponents of spacetime tensors on a particular coordinate 
system; they thus transform under spacetime coordinate 
transformations. They run over 0,1,2,3. Lower case 
Greek indices a, /?, 7 . . . from the start of the alphabet la- 
bel position or momentum coordinates on 8 dimensional 
phase space that are not associated with coordinates on 
spacetime. They run over 0, 1,2,3 and do not transform 
under spacetime coordinate transformations. In Sec. fVl 
and just in that section, indices a, /?, 7, (5, e, . . . from the 
start of the Greek alphabet run over 1 . . . iV, and indices 
A, fi, V, p,a, . . . from the second half of the alphabet run 
over 1 . . .M . Bold face quantities generally denote vec- 
tors, as in J = (Ji, . . . , Jm), although in SecIUthe bold 
faced notation is used for differential forms. 



^ In coordinates t = t — r, r, 9, (p, the explicit form of the asym 
totic solution can be obtained by taking Eq. (3.1) of Ref. 
eliminating the phases Xlmkn using Eq. (8.29) of Ref. f68|, and 
making the identifications qr = Qr[t — r — ta + ir{—\-o) — 
*9(--^6lo)] - TTrArO, qo = f^6) [< " ^ -Jo + ir(-Xro) - tg {- Xgo)] - 

TeXeo, and = n^[t - r - to + tr(-\ro) - ki-Xgo)] + <t>o - 
(f>r{-Xro) + 4>e{-Xeo)- 

The function T^p depends on q^ and 4> only through the combi- 
nation q^ — (f>. This allows us to show that the two-timescale form 
l|1.8|l of the metric reduces to a standard Taylor series expansion, 
locally in time near almost every value to of i. For equatorial 
orbits there is no dependence on qg, and the e dependence of 
the metric has the standard form up to linear order, in coordi- 
nates {t',r',e',(j>') defined by t' = {t - io)/e + [fi^°\io)/£]/i^rO, 
4,' = 4, + u>^o[fr''\iQ)/e\/u>rQ - [ff\io) I s], r' = r, e' = 9, 

where Wr-o = fr^^' {to), lj^o = f^^' {io)^ and for any number x, 
[x\ = x + 2iTn where the integer n is chosen so that < [a;] < 27r. 
A similar construction works for circular orbits for which there is 
no dependence on qr- For generic orbits a slightly more involved 
construction works, but only if uJro/<^4>o is irrational [8d |. which 
occurs for almost every value of to. 



II. EXTREME-MASS RATIO INSPIRALS IN 
KERR FORMULATED USING ACTION-ANGLE 
VARIABLES 

In this section we derive the form of the fundamental 
equations describing the inspiral of a point particle into 
a Kerr black hole, using action-angle type variables. Our 
final result is given in Eqs. (|2.47p below, and the prop- 
erties of the solutions of these equations are analyzed in 
detail in the remaining sections of this paper. 

The description of geodesic motion in Kerr in terms 



93|, 
94|. 



of action angle variables was first given by Schmidt 
and has been reviewed by Glampedakis and Babak 
We follow closely Schmidt's treatment, except that we 
work in an eight dimensional phase space instead of a 
six dimensional phase space, thus treating the time and 
spatial variables on an equal footing. We also clarify the 
extent to which the fundamental frequencies of geodesic 
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motion are uniquely determined and gauge invariant, as 
claimed by Schmidt. 

We start in subsection III Al by reviewing the geometric 
definition of action angle variables in Hamiltonian me- 
chanics, which is based on the Liouville- Arnold theorem 
[9^. This definition docs not apply to geodesic motion 
in Kerr, since the level surfaces defined by the conserved 
quantities in the eight dimensional phase space are non- 
compact. In subsection III Bl we discuss how generalized 
action angle variables can be defined for non-compact 
level surfaces, and in subsection III CI we apply this to 
give a coordinate-independent construction of general- 
ized action angle variables for generic bound geodesies in 
Kerr. Subsection III Dl specializes to Boyer-Lindquist co- 
ordinates on phase space, and describes explicitly, follow- 
ing Schmidt [93|, the explicit canonical transformation 
from those coordinates to the generalized action angle 
variables. 

We then turn to using these variables to describe a 
radiation-reaction driven inspiral. In subsection III El we 
derive the equations of motion in terms of the general- 
ized action angle variables. These equations define a flow 
on the eight dimensional phase space, and do not explic- 
itly exhibit the conservation of rest mass. In subsection 
IIIFI we therefore switch to a modified set of variables 
and equations in which the conservation of rest mass is 
explicit. We also augment the equations to describe the 
backreaction of gravitational radiation passing through 
the horizon of the black hole. 



A. Review of action-angle variables in geometric 
Hamiltonian mechanics 

We start by recalling the standard geometric frame- 
work for Hamiltonian mechanics [osj . A Hamiltonian 
system consists of a 2A^-dimensional differentiable man- 
ifold A4 on which there is defined a smooth function 
H (the Hamiltonian), and a non-degenerate 2-form Qab 
which is closed, W[a^bc] — 0. Defining the tensor fl'^^ by 
il'^''flbc = S^, the Hamiltonian vector field is defined as 

v" = f7"''VbiJ, (2.1) 

and the integral curves of this vector fields give the mo- 
tion of the system. The two form fiab is called the sym- 
plectic structure. Coordinates {qa,Pa) with 1 < a < 
are called symplectic coordinates if the symplectic struc- 
ture can be written as = dpa A dqa, i.e. flab — 

We shall be interested in systems that possess A — 1 
first integrals of motion which, together with the Hamil- 
tonian H, form a complete set of A^ independent first 
integrals. We denote these first integrals by Pq, 1 < a < 
N, where Pi = H . These quantities are functions on 
for which the Poisson brackets 

{P^,H} = n''\VaPc.){VbH) (2.2) 



vanish for 1 < a < A. If the first integrals satisfy the 
stronger condition that all the Poisson brackets vanish, 

{P„,P^} = (2.3) 

for 1 < a, /3 < A^, then the first integrals are said to 
be in involution. If the 1-forms V aPa for 1 < a < A" 
are linearly independent, then the first integrals are said 
to be independent. A system is said to be completely 
integrable in some open region U in M if there exist A^ 
first integrals which are independent and in involution at 
every point of lA. 

For completely integrable systems, the phase space A4 
is foliated by invariant level sets of the first integrals. For 
a given set of values p = {pi, . . . ,pn), we define the level 
set 

Mp = {xeM\Paix)^Pa,l<a<N}, (2.4) 

which is an A^-dimensional submanifold of A4 . The level 
sets are invariant under the Hamiltonian fiow by Eq. 
(12. 2p . Also the pull back of the symplectic structure ft 
to A^p vanishes, since the vector fields Va defined by 

< = n'^'^^bPa (2.5) 

for 1 < a < A^ form a basis of the tangent space to A^p 
at each point, and satisfy r^afow" wj^ = for 1 < a, /3 < A 
by Eq. (lOll . 

A classic theorem of mechanics, the Liouville- Arnold 
theorem [95|, applies to systems which are completely 
integrable in a neighborhood of some level set A4p that 
is connected and compact. The theorem says that 

• The level set A^p is diffeomorphic to an A-torus 
r^. Moreover there is a neighborhood V of A^p 
which is diffeomorphic to the product x B where 
B is an open ball, such that the level sets are the 
A^-tori. 

• There exist symplectic coordinates {qa, J a) for 1 < 
a < N (action-angle variables) on V for which the 
angle variables qa are periodic, 

qa+2TT= qa, 

and for which the first integrals depend only on the 
action variables, Pa = Pq(Ji, . . . , Jn) for 1 < a < 
N. 

An explicit and coordinate-invariant prescription for 
computing a set of action variables Ja is as follows 95] . 
A symplectic potential is a 1-form which satisfies 
d& = n. Since the 2-form Jl is closed, such 1-forms 
always exist locally. For example, in any local symplec- 
tic coordinate system {qa,Pa), the 1-form — Padqa is 
a symplectic potential. It follows from the hypotheses of 
the Liouville-Arnold theorem that there exist symplec- 
tic potentials that are defined on a neighborhood of A^p 
[9^. The first homotopy group Hi(A^p) is defined to 
be the set of equivalence classes of loops on A^p, where 
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two loops are equivalent if one can be continuously de- 
formed into the other. Since A4p is diffeomorphic to the 
A''-torus, this group is isomorphic to (Z^, +), the group 
of iV-tuples of integers under addition. Pick a set of 
generators 71, . . . ,7Ar of ni(A^p), and for each loop 7^ 
define 



1 

ZTT 



0. 



(2.6) 



This integral is independent of the choice of symplectic 
potential 0.^^ It is also independent of the choice of loop 
-fa in the equivalence class of the generator of Hi (Alp), 
since if 7q and 7^ are two equivalent loops, we have 



/ = 



on 



n 



& = d& = fl = 0. (2.7) 



n 



Here 7?, is a 2-dimensional surface in A4p whose boundary 
is la — I'a^ we have used Stokes theorem, and in the last 
equality we have used the fact that the pull back of $7 to 
the level set A^p vanishes. 

Action-angle variables for a given system are not 
unique [13] ■ There is a freedom to redefine the coor- 
dinates via 



qa 



(2.8) 



where Aaf3 is a constant matrix of integers with deter- 
minant ±1, and AapBa-y — Spy. This is just the free- 
dom present in choosing a set of generators of the group 
ni(Alp) - (Z^,-H). Fixing this freedom requires the 
specification of some additional information, such as a 
choice of coordinates on the torus; once the coordinates 
Qa are chosen, one can take the loops to be the curves 
qi3 = constant for /3 ^ a. There is also a freedom to re- 
define the origin of the angle variables separately on each 
torus: 



qa 



dZjJp) 

dJa 



Jot ^ Ja- 



(2.9) 



Here Z[Jp) can be an arbitrary function of the action 
variables. 



B. Generalized action-angle variables for 
non-compact level sets 

One of the crucial assumptions in the Liouville- Arnold 
theorem is that the level set Alp is compact. Unfortu- 
nately, this assumption is not satisfied by the dynami- 
cal system of bound orbits in Kerr which we discuss in 
Sec. IH CI below, because we will work in the 8 dimen- 
sional phase space and the motion is not bounded in the 



The type of argument used in Ref. 1991 can be used to show that 
the puUback to A^p of the difference between two symplectic 
potentials is exact since it is closed. 



time direction. We shall therefore use instead a general- 
ization of the Liouville- Arnold theorem to non-compact 
level sets, due to Fiorani, Giachetta and Sardanashvily 

Consider a Hamiltonian system which is completely 
integrable in a neighborhood U & connected level set 
A^p, for which the N vector fields (|2.5p are complete on 
U, and for which the level sets A^p' foliating U are all 
diffeon iorp hic to one another. For such systems Fiorani 
et. al. [96| prove that 

• There is an integer k with < k < N such that 
the level set A^p is diffeomorphic to the product 
T'^ X R^"*^, where R is the set of real numbers. 
Moreover there is a neighborhood V of A^p which 
is diffeomorphic to the product T'"' x R^^*"' x B 
where B is an open ball. 

• There exist symplectic coordinates {qa, Ja) for 1 < 
a < N (generalized action-angle variables) on V for 
which the first k variables qa are periodic. 



2n : 



qa, 



l<a<k, 



and for which the first integrals depend only on the 
action variables, Pa — Pa{Ji, ■ ■ ■ , Jn) for 1 < a < 
N. 

Thus, there are k compact dimensions in the level sets, 
and N^k non-compact dimensions. In our application to 
Kerr below, the values of these parameters will be fc = 3 
and N -k = l. 

The freedom in choosing generalized action-angle vari- 
ables is larger than the corresponding freedom for action- 
angle variables discussed above. The first k action vari- 
ables can be computed in the same way as before, via the 
integral (|2.6p evaluated on a set of generators 71 , . . . , 7fc 
of Hi(A4p), which in this case is isomorphic to (Z*^,-!-). 
This prescription is unique up to a group of redefinitions 
of the form [cf. Eq. above] 



k 

E 

/3=1 



Aal3qf3, Ja ^'^BapJ/S, (2-10) 

/3=1 



for 1 < a < fc, where the k x k matrix Aaf3 is a 
constant matrix of integers with determinant ±1, and 
AapBa-y — Spy. There is additional freedom present in 
the choice of the rest of the action variables Jk+i, ■ ■ ■ , Jn- 
As a consequence, the remaining freedom in choosing 
generalized action-angle variables consists of the trans- 
formations (12. 9p discussed earlier, together with trans- 
formations of the form 



qa 



Aapqp, Ja — > BapJp, 



(2.11) 



where Aap and Bap are constant real N x N matrices 
with AapBa-f — Spy such that Ji, . . . ,Jk are preserved. 

In generalized action-angle variables, the equations of 
motion take the simple form 



dH(3) 

dJa 



(2.12) 
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and 



We define the quantities 



dH{3) 

dJa 



(2.13) 



(2.14) 



which are angular frequencies for 1 < < fc but not for 
k + 1 < a < N. The solutions of the equations of motion 
are then 



qa(t) = Cla{3o)t + qaO 
Ja{t) — JaO, 



(2.15a) 
(2.15b) 



for some constants Jo and qo. 



C. Application to bound geodesic motion in Kerr 

We now apply the general theory discussed above to 
give a coordinate-invariant definition of action-angle vari- 
ables for a particle on a bound orbit in the Kerr space- 
time. We denote by {MK,gab) the Kerr spacetime, and 
we denote by and rj"" the timelike and axial Killing 
vector fields. The cotangent bundle over A^k forms an 
8-dimensional phase space Ai = T*A4k- Given any co- 
ordinate system x'^ on the Kerr spacetime, we can define 
a coordinate system {x'^,pj,) on M, such that the point 
{x'^,p^) corresponds to the covector or one form p^dx'^ 
at cc'^ in A^K- The natural symplectic structure on 
is then defined by demanding that all such coordinate 
systems {x'^ ^p^) be symplectic [95|. The Killing vector 
fields and ry" on A^k have natural extensions to vec- 
tor fields on phase space which Lie derive the symplectic 
structure. 

Consider now a particle of mass /i on a bound geodesic 
orbit. A Hamiltonian on that generates geodesic mo- 
tion is given by 

H{x'',Pu) = \9''%xP)p,p.- (2.16) 

this definition is independent of the choice of coordinate 
system x'^ . If we interpret Pi, to be the 4- momentum of 
the particle, then the conserved value of H is — /i^/2, and 
the evolution parameter is the afhne parameter X — t/ fi 
where t is proper time. 

As is well known, geodesies on Kerr possess three first 
integrals, the energy E = — CPa, the z-component of 
angular momentum — rj'^Pa , and Carter constant Q — 
Q'^^PaPb where Q"'' is a Killing tensor 98] . Together with 
the Hamiltonian we therefore have four first integrals: 



P„ = {Po,Pi,P2,P3) = {H,E,L,,Q). 



(2.17) 



An explicit computation of the 4-form dH f\dE f\dL ^ AdQ 
on M shows that it is non vanishing for bound orbits ex- 
cept for the degenerate cases of circular (i.e. constant 



Boyer-Lindquist radial coordinate) and equatorial or- 
bits. Also the various Poisson brackets {Pa,Pi3} vanish: 
{E,H} and {Lz;H} vanish since and ry" are Killing 
fields, {E, Lz} vanishes since these Killing fields com- 
mute, {Q, H} vanishes since is a Killing tensor, and 
finally {i?, Q} and {Lz,Q} vanish since the Killing ten- 
sor is invariant under the flows generated by and 77". 
Therefore for generic orbits the theorem due to Fiorani 
et. al. discussed in the last subsection applies. The rel- 
evant parameter values are k = 3 and TV = 4, since the 
level sets Mp are non-compact in the time direction only. 
Thus geodesic motion can be parameterized in terms of 
generalized action-angle variables. 

We next discuss how to resolve in this context the 
non-uniqueness in the choice of generalized action an- 
gle variables discussed in the last subsection. Consider 
first the freedom (|2.10|) associated with the choice of gen- 
erators of ni(A^p). One of these generators can be cho- 
sen to be an integral curve of the extension to A4 of 
the axial Killing field 77". The other two can be chosen 
as follows. Let tt : AI — > A4k be the natural projec- 
tion from phase space Ai to spacetime AIk that takes 
(x^^pij) to x'^. A loop (a;''(A),p^(A)) in the level set Mp 
then projects to the curve x'^{X) in n(A4p). Requiring 
that this curve intersect the boundary of 7r(Alp) only 
twice determines the two other generators of Hi (Alp). ^■^ 
The resulting three generators coincide with the gen- 
erators obtained from the motions in the r, 9 and (jj 
directions in Boyer-Lindquist coordinates [93]. We de- 
note the resulting generalized action-angle variables by 

The remaining ambiguity (|2.1ip is of the form 



(2.18) 



where i runs over r, and 4> and the parameters 7 and 
are arbitrary. The corresponding transformation of the 
frequencies (|2.14p is 



7 



7 



(2.19) 



A portion of this ambiguity (the portion given by 7 = 1, 
v"^ = ^ 0) is that associated with the choice of rota- 
tional frame, (j) ^ (j) + U,t where £7 is an angular velocity. 
It is not possible to eliminate this rotational-frame am- 
biguity using only the spacetime geometry in a neighbor- 
hood of the orbit. In this sense, the action angle variables 
are not uniquely determined by local geometric informa- 
tion. However, we can resolve the ambiguity using global 
geometric information, by choosing 



Jt - — / 

27r 



(2.20) 



It 



One can check that the two other assumptions in the theorem 
listed in the second paragraph of Sec. Ill Bl are satisfied. 
This excludes, for example, loops which wind around twice in 
the r direction and once in the 6 direction. 
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where 74 is an integral curve of length 27r of the extension 
to M of the timehke KiUing field C^.^'' The definition 
()2.20p is independent of the choice of such a curve -ft and 
of the choice of symplectic potential 0. 

To summarize, we have a given a coordinate- 
invariant definition of the generalized action-angle vari- 
ables {qt, qr,qe, q(t,, Jt, Jr, J9, J(t>) for generic bound orbits 
in Kerr. These variables are uniquely determined up to 
relabeling and up to the residual ambiguity (|2.9I1 . A 
similar construction has been given by Schmidt [93;], ex- 
cept that Schmidt first projects out the time direction 
of the level sets, and then defines three action variables 
{Jr, Jg, J4,) and three angle variables (g^i <le, <l4>)- 



D. Explicit expressions in terms of 
Boyer-Lindquist coordinates 

In Boyer-Lindquist coordinates {t^r,9,(j)), the Kerr 
metric is 



2Mr\ , o S 



dr"^ + E de'^ 



2Ma^r 



r + a 
AMar 



SIT? 6 1 sin^ I 



■ sin^ 6 dt 



where 

S = r- 



2,2 2/1 

' a cos 0, 



(2.21) 



r'^-2Mr + a^, (2.22) 



and M and a are the black hole mass and spin parame- 
ters. The timelike and axial Killing fields are = d/dt 
and ff — d/d(f), and so the energy and angular momentum 
are 



E 



-Pt 



and 



The Carter constant is given by [o^ 
Q = pI + a^cos^e* (//^ -p?) 
and the Hamiltonian (|2.16p is 



(2.23a) 



(2.23b) 



cot^Spi (2.23c) 



H 



a sm 



-Pr 



2E"'' ■ 21]"" ■ 2^sm^0 
1 2 



2EA 



(2.23d) 



Following Schmidt [9^, we can obtain an invertible 
transformation from the Boyer-Lindquist phase space co- 
ordinates {x'^ ,pv) to the generalized action angle vari- 
ables {qa,Ja) as follows. Equations (|2.23p can be in- 
verted to express the momenta Pi, in terms of x"^ and the 
four first integrals 



1 9 
--n\E,Lz,i 



(2.24) 



up to some signs [98|: 



Pt = -E, P4, = L2, Pr = ±— pe = ±\/V0{e). 



A 

Here the potentials Vr{r) and Ve{9) are defined by 

1 2 



(2.25) 



Vr{r) = [{r^ + a'^)E ~ aL, 

-A [fi^r^ + {L, - aEf + Q] , (2.26a) 



Ve{e) ^ Q- 



ili' - E')a' + 



sm 



cos^6'.(2.26b) 



Using these formulae together with the symplectic poten- 
tial — Pi^dx" in the definitions (|2.6p and (|2.20p gives 





2^. 


(f — r—dr 
T A 


(2.27a) 


Je = 




<j> VVede 


(2.27b) 


J4, — 


2^. 


j> p^dcj) = 


(2.27c) 




2^. 






Jt = 


/ ptdt = -E. 
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(2.27d) 



These expressions give the action variables as functions of 
the first integrals, Jq, — Ja{Pp)- The theorem discussed 
in Sec. Ill Bl above guarantees that these relations can be 
inverted to give 



Pa (J/3). 



(2.28) 



Next, to obtain expressions for the corresponding gen- 
eralized angle variables, we use the canonical transforma- 
tion from the symplectic coordinates {x'^ ,pi,) to {qa, J a) 
associated with a general solution of the Hamilton Jacobi 
equation 



H 



dS 

dx'' 



dS 



(2.29) 



As shown by Carter 98], this equation is separable and 
the general solution^ ^ can be written in terms of the first 



The Killing field encodes global geometric information since 
it is defined to be timelike and of unit norm at spatial infinity. 



^ As indicated by the it signs in Eq. I|2.3ip . there are actually four 
different solutions, one on each of the four coordinate patches 
on which {x'^ , Pa) are good coordinates, namely sgn(pr) = il, 
sgn(pe) = ±1. 
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integrals Pa 

S{x'',Pa,X) = -HX + W{x\Pa) (2.30) 
where H = — /i^/2, 

yV{x\Pa) = -Et + L^cj) ± Wr{r) ± Weie), (2.31) 



A ' 



and 



Wr{r) = / dr 



We{0) = J d9^g. 



(2.32) 



(2.33) 



Using the relation (|2.28p the function W can be expressed 
in terms of the Boyer-Lindquist coordinates x'^ and the 
action variables Ja- 



W^W{x\Ja). 



(2.34) 



This is a type II generating function that generates 
the required canonical transformation from {x'^,p^) to 

{Qa ^ Jol) '- 



Pv 



dx" 
dW 

dJa 



(x^J^) 
(x^J^). 



(2.35a) 
(2.35b) 



Equation (|2.35ap is already satisfied by virtue of the def- 
inition ((2?3T|l of W together with Eqs. ((2?25)) . Equation 
(j2.35b[) furnishes the required formulae for the general- 
ized angle variables Qa-^^ 

Although it is possible in principle to express the first 
integrals Pa in terms of the action variables Ja using Eqs. 
(|2.27p , it is not possible to obtain explicit analytic expres- 
sions for Pa{J())- However, as pointed out by Schmidt 
[9^, it is possible to obtain explicit expressions for the 
partial derivatives dPa/dJp, and this is sufhcient to com- 
pute the frequencies Via- We review this in appendix 1X1 



E. Application to slow inspiral motion in Kerr 

The geodesic equations of motion in terms of the gen- 
eralized action angle variables (ga, Ja) are [cf. Eqs. (|2.12p 
- ((TTi)) above] 



dqa 
dX 

dJa 

ItX 



0, 



(2.36a) 
(2.36b) 



The freedom 1)2. 9p to redefine the origin of the angle variables on 
each torus is just the freedom to add to W any function of Pa. 
We choose to resolve this freedom by demanding that = at 
the minimum value of r, and = at the minimum value of 9. 



for < a < 3. Here X = r/fi where t is proper time and 
fi is the mass of the particle. In this section we derive the 
modifications to these equations required to describe the 
radiation-reaction driven inspiral of a particle in Kerr. 
Our result is of the form 

^ = naiJf3)+tl\fa{qf3,Jf3), (2.37a) 

^ = ti^Fa{qp,Jp). (2.37b) 

We will derive explicit expressions for the forcing terms 
fa and Fa in these equations. 

The equation of motion for a particle subject to a self- 
acceleration a" is 



dA2 



dx" dxP 



2 1/ 

/I a . 



(2.38) 



Rewriting this second order equation as two first order 
equations allows us to use the Jacobian of the coordinate 
transformation {x'^,p^} — > {qa, Ja} to relate the forcing 
terms for the two sets of variables: 



dx" 
H 
dpi, 
dX 



9"'' P., 



(2.39a) 
(2.39b) 



We start by deriving the equation of motion for the 
action variables Ja- Taking a derivative with respect to 
A of the relation Ja = Ja{x^ tPv) and using Eqs. (|2.39|) 
gives 



dJa 



P 



dJg 
dx" 

dJg 

dx" 

B T 
+H - — Oj, 

Op,y 



9 Pa 



dJa dpi, 
dpi, dX 
1 dJa 
28^ 



9''%P<yPp 



(2.40) 



The term in square brackets must vanish identically since 
J a is conserved in the absence of any acceleration a^. 
Rewriting the second term using J a — JaiPp) and the 
chain rule gives an equation of motion of the form (j2.37b[) . 
where the forcing terms Fa are 



Fa 



dJa ( dPp 
dPp \ dp. 



(2.41) 



Here the subscript x on the round brackets means that 
the derivative is to be taken holding x'^ fixed. When the 
sum over (3 is evaluated the contribution from Pi = H 
vanishes since a^p" = 0, and we obtain using Eqs. (|2.17p 
and ^(m\ 



Ft 

Fr 

Fe 
F^ 



at, 



dJr 

'dE 
dJe_ 

' dE 



at 



at 



dJr 
dJo 



dQ 



aq 



dJr 
dL 
dJe 



-acjy, 



(2.42a) 
(2.42b) 

(2.42c) 
(2.42d) 
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Here we have defined oq = 2Q"'^p^aa and the various 
coefBcients dJa/dPp are given expHcitly as functions of 
Pa in Appendix lAl 

We use a similar procedure to obtain the equation of 
motion (|2.37ap for the generahzed angle variables qa ■ Dif- 
ferentiating the relation q^. = qa{x'^ iPv) with respect to 
A and combining with the two first order equations of 
motion (|2.39p gives 



dqg 
dX 



dqa 

i^qa 

-[I - — a^. 



(2.43) 



By comparing with Eq. (|2.36ap in the case of vanishing 
acceleration we see that the term in square brackets is 
^aiJp)- This gives an equation of motion of the form 
(|2.37ap . where the where the forcing term /q, is 



Ja = — ay. 



(2.44) 



Using the expression (|2.35b|l for the angle variable qo 
together with Jq, — JaiPji) gives 



dqg 
dpi, 



dPy 
fdW\ 



'dPp ( 




_dJa \ 




d 


(dp,y\ 


.9P, 


\dJcJ\ 



(2.45) 



This yields for the forcing term 
'dP^ 



fa 



7 

dpi, 



m 

dJa 



(. 



ydPsdP^ 

2 



(dW\ dPp d^Je 
\dPp),^ dJ, dP-,dPs 



(2.46) 



In this expression the first two factors are the same as 
the factors which appeared in the forcing term (|2.4ip for 
the action variables. The quantities dPs/dJa, dPp/dJe 
and d'^Je/{dPjdPs) can be evaluated explicitly as func- 
tions of Pa using the techniques discussed in Appendix 
El The remaining factors in Eq. (|2.46p can be evaluated 
by differentiating the formula (|2.3ip for Hamilton's prin- 
cipal function W and using the formulae (|2.26p for the 
potentials Vr and Vg. 



Rescaled variables and incorporation of 
backreaction on the black hole 



We now augment the action-angle equations of motion 
p.37p in order to describe the backreaction of the gravi- 
tational radiation on the black hole. We also modify the 
equations to simplify and make explicit the dependence 
on the mass /i of the particle. The resulting modified 



equations of motion, whose solutions we will analyze in 
the remainder of the paper, are 



dqg 
dT 



= uja{Pj,MB) + eg^^\qA, Pj,Mb) 

+e^9^^HqA,Pr.MB) + 0(e3), (2.47a) 



^ = eG\'\qA,P„MA) + e'G^;'\qA,P„MB) 



dMA 
dT 



+0{e^), (2.47b) 
e^GAiqA,Pj,MB) + Oie''). (2.47c) 



Here a runs over 0,1,2,3, i, j run over 1,2,3, A, B 
run over 1,2, qA = {qr,qs), Ma = (Mi,M2) and Pi = 

(Pi, P2, Ps)- Also all of the functions w^, ga\ ga\ G^^\ 
(2'] 

Gi and Ga that appear on the right hand sides are 
smooth functions of their arguments whose precise form 
will not be needed for this paper (and are currently un- 
known aside from uja)- 

Our final equations (|2.47p are similar in structure to 
the original equations (j2.37p . but there are a number of 
differences: 

• We have switched the independent variable in the 
differential equations from affinc parameter A to 
proper time t = //A. 

• We have introduced the ratio 



M 



of the particle mass and black hole mass M , and 
have expanded the forcing terms as a power series 
in e. 

• The forcing terms g^a \ ga \ Gf'\ g[^\ and Ga de- 
pend only on the two angle variables qA = (9r, 9e)j 
and are independent of qt and q^. 

• Rather than evolving the action variables Ja, we 
evolve two different sets of variables. Pi and Ma- 
The first of these sets consists of three of the first 
integrals of the motion, with the dependence on the 
mass /i of the particle scaled out: 

A - {Pi,P2,Pz) = {E/^l,LJ^l,Q/^?). (2.49) 

The second set consists of the mass and spin pa- 
rameters of the black hole, which gradually evolve 
due to absorption of gravitational radiation by the 
black hole: 



Ma = (Mi,M2) = (Af,a). 



(2.50) 



We now turn to a derivation of the modified equations 
of motion (|2.47p . The derivation consists of several steps. 
First, since the mapping (|2.27p between the first integrals 
Pa and the action variables Ja is a bijection, we can use 
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the Pa as dependent variables instead of Jq.^^ Equation 
(|2.37ap is unmodified except that the right hand side is 
expressed as a function of Pa instead of Ja- Equation 
(|2.37bp is replaced by 



dPa 



dPg 



(2.51a) 



Second, we switch to using modified versions Pa of 
the first integrals Pa with the dependence on the mass fi 
scaled out. These rescaled first integrals are defined by 



Pa = {H,E,L,,Q) 

= [H|^x\E/^i,L,/^Ji,Q/^?). 



(2.52) 



We also change the independent variable from afhnc pa- 
rameter A to proper time r = /iA. This gives from Eqs. 
(|2.37p and (|2.44p the system of equations 



dPg 

dT 



dPg 



(2.53a) 
(2.53b) 



where we have defined — (2, 1, 1, 2). 

Third, we analyze the dependence on the mass fi of 
the right hand sides of these equations. Under the trans- 
formation (x'^jPi,) {x'^ , sp^) for s > 0, we obtain 
the following transformation laws for the first integrals 
(|2.24p . the action variables (|2.27p . and Hamilton's prin- 
cipal function (|2.3ip : 

Pa s''"Pa with = (2,1,1,2), (2.54a) 
Ja sJa, (2.54b) 
W -> sW. (2.54c) 



From the definitions (|2.14p and (j2.35bp of the angular 
frequencies Q.a and the angle variables qa we also deduce 



Via 



silo 

qa- 



(2.55a) 
(2.55b) 



Similarly, if we write the angle variable qa as a func- 
tion qa{x'^ ,Piy) of x^ and pi,, then the scaling law (|2.55b[) 
implies that qa{x^ , spu) — qa{x'^,Pv), and it follows that 
the coefficient of the 4-acceleration in Eq. (|2.53ap is 



dqa, cr N dqa, „ , dqa 
opy opi, opi, 



(a;",u,), (2.57) 



where is the 4-velocity. This quantity is also indepen- 
dent of /i at fixed P/s. We will denote this quantity by 
failp^Pp)- It can be obtained explicitly by evaluating 
the coefficient of ai, in Eq. (|2.46p at Pa — Pa, Pv = Uy. 
A similar analysis shows that the driving term on the 
right hand side of Eq. (j2.53bp can be written in the form 



The resulting rescaled equations of motion are 



dqa 
dT 



t^aiPp) + fa{<lti,Pl3)au, 



^ = Fa{<l!}^P0)au- 



(2.58) 

(2.59a) 
(2.59b) 



Note that this formulation of the equations is completely 
independent of the mass /i of the particle (except for the 
dependence on of the radiation reaction acceleration 
ai, which we will discuss below). 

Fourth, since Pq = H = — /i^/2, the rescaled variable 
is Po = -1/2 from Eq. ([232]) . Thus we can drop the 
evolution equation for Pq, and retain only the equations 
for the remaining rescaled first integrals 



P, = (Pi,P2,P3) - iE,L,,Q). 



(2.60) 



We can also omit the dependence on Pq in the right hand 
sides of the evolution equations (|2.59p , since Pq is a con- 
stant. This yields 



dqa 

dT 



UJa{Pj)+ fa{<lf},Pj)au, 



(2.61a) 
(2.61b) 



If we write the angular velocity JIq, as a function uja{Pp) 
of the first integrals P^, then it follows from the scalings 
(I2.54ap and (|2.55ap that the first term on the right hand 
side of Eq. (|2.53ap is 



na _ UJa{Pp) _ tOaill^^Pp) 



= UJa{Pl3)- (2.56) 



This quantity is thus independent of /i at fixed P^, as we 
would expect. 



Fifth, the self-acceleration of the particle can be ex- 
panded in powers of the mass ratio s = fJ./M as 



+ e2a(2)+0(£3). 



(2.62) 



Here a[}^ is the leading order self-acceleration derived by 
Mino, Sasaki and Tanaka [3| and by Quinn and Wald 
[45l |. discussed in the introduction. The subleading self- 
acceleration a^^ has been computed in Refs. [H, [H, [sgI . 



Note that fid/dpu cannot be simplified to d/dui, because we are 
^'^ Note that since the variables J a are adiabatic invariants, so are working in the eight dimensional phase space M where /i is a 

the variables Pa- coordinate and not a constant. 
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[13, nil- The accelerations al}^ and al?-* are independent 
of fi and thus depend only on and u^, or, equivalently, 
on g„ and Pi. This yields the system of equations 



-0(e3) 



Here the forcing terms are given by 



G 



is) 



(2.63a) 

(2.63b) 

(2.64a) 
(2.64b) 



for s = 1,2. 

The formula (|2.62p for the self-acceleration, with the 

explicit formula for al^^ from Refs. [3|4^, is valid when 
one chooses the Lorentz gauge for the metric pertur- 
bation. The form of Eq. (j2.62l) is also valid in a vari- 
ety of other gauges; see Ref. [99] for a discussion of the 
gauge transformation properties of the self force. How- 
ever, there exist gauge choices which are incompatible 
with Eq. ([2^ . which can be obtained by making e- 
dependent gauge transformations. We shall restrict at- 
tention to classes of gauges which are consistent with our 
ansatz (jl.Sp for the metric, as discussed in Sec. II El above. 
This class of gauges has the properties that (i) the devia- 
tion of the metric from Kerr is < e over the entire inspi- 
ral, and (ii) the expansion (|2.62p of the self-acceleration is 
valid. These restrictions exclude, for example, the gauge 
choice which makes al}^ = 0, since in that gauge the 
particle does not inspiral, and the metric perturbation 
must therefore become of order unity over an inspiral 
time. We note that alternative classes ofgauees have 
been suggested and explored by Mino [s^, [72, [Sll, [13 ■ 

Sixth, from the formula (j2.35bp for the generalized an- 
gle variables qa together with Eqs. (|2.3ip and (|2.27d[) it 
follows that qt can be written as 



qt=t + ft{r,e,P^) 



(2.65) 



for some function ft. All of the other angle and action 
variables are independent of t. Therefore the vector field 
d/dt on phase space is just d/dqt, the symmetry t — > 
t -(- At with x% fixed is the same as the symmetry 
qt ^ qt + with qr,qe, q<p and Ja fixed. Since the self- 
acceleration as well as the background geodesic motion 
respect this symmetry, all of the terms on the right hand 
side of Eqs. (|2.63p must be independent of g*. A similar 
argument shows that they are independent of q^. This 
gives 



dqa 
dr 



^ = Lu^iP,) + egi'\qA,P,)+e^g^^\qA,P,) 



^ = eGf\^qA.P,)+e'Gf\qA,P,) 



(2.66a) 
(2.66b) 



where qA = {qr,qd)- 

Seventh, consider the evolution of the black hole back- 
ground. So far in our analysis we have assumed that the 
particle moves in a fixed Kerr background, and is subject 

to a self-force = eai^' -I- e'^al?^ + 0{e'^). In reality, the 
center of mass, 4-momentum and spin angular momen- 
tum of the black hole will gradually evolve due to the 
gravitational radiation passing through the event hori- 
zon. The total change in the mass M of the black hole 
over the inspiral timescale ~ M/e is ~ Me. It follows 
that the timescale for the black hole mass to change by 
a factor of order unity is ~ Mje^ . The same timescale 
governs the evolution of the other black hole parameters. 

This effect of evolution of the black hole background 
will alter the inspiral at the first subleading order (post-1- 
adiabatic order) in our two-timescale expansion. A com- 
plete calculation of the inspiral to this order requires solv- 
ing simultaneously for the motion of the particle and the 
gradual evolution of the background. We introduce the 
extra variables 



MA = (Mi,M2) = (M,a), 



(2.67) 



the mass and spin parameters of the black hole. We 
modify the equations of motion (j2.66p by showing explic- 
itly the dependence of the frequencies uja and the forcing 
functions g^a^ and g["^ on these parameters (the depen- 
dence has up to now been implicit). We also add to the 
system of equations the following evolution equations for 
the black hole parameters: 



dMA 

dr 



^GA{qB,Pj.MB) + 0[e^), 



(2.68) 



where A — 1,2. Here Ga are some functions describing 
the fluxes of energy and angular momentum down the 
horizon, whose explicit form will not be important for 
our analyses. They can in principle be computed using, 
for example, the techniques developed in Ref. 100].^^ 
The reason for the prefactor of is that the evolution 
timescale for the black hole parameters is ~ Mje^, as 
discussed above. The functions G a are independent of 
qt and q^ for the reason discussed near Eq. (|2.66p : the 
fluxes through the horizon respect the symmetries of the 
background spacetime. Finally, we have omitted in the 
set of new variables (|2.67p the orientation of the total 



These techniques naturally furnish the derivatives of Mj\ with 
respect to Boyer Lindquist time t, not proper time r as in Eq. 
II2.()8| I. However this difference is unimportant; one can easily 
convert from one variable to the other by multiplying the func- 
tions Ga by the standard expression for dt/dr 101], 



dt 



where = \^r^~+~a^. This expression can be written in terms 
of of QA , Pi and Ma , and is valid for accelerated motion as well 
as geodesic motion by Eqs. 112.251 1 and I l2.39al l. 
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angular momentum, the location of the center of mass, 
and the total linear momentum of the system, since these 
parameters are not coupled to the inspiral motion at the 
leading order. However, it would be possible to enlarge 
the set of variables Ma to include these parameters with- 
out modifying in any way the analyses in the rest of this 
paper. 

These modifications result in the final system of equa- 
tions jMll). 

Finally we note that an additional effect arises due to 
the fact that the action-angle variables we use are de- 
fined, at each instant, to be the action-angle variables 
associated with the black hole background at that time. 
In other words the coordinate transformation on phase 
space from (a;'^,py) —^ {qa, Jo) acquires an additional de- 
pendence on time. Therefore the Jacobian of this trans- 
formation, which was used in deriving the evolution equa- 
tions (j2.37p . has an extra term. However, the correspond- 
ing correction to the evolution equations can be absorbed 

(2) 

into a redefinition of the forcing term ga ■ 

G. Conservative and dissipative pieces of the 
forcing terms 

In this subsection we define a splitting of the forc- 
ing terms Qa and Gi in the equations of motion (|2.47p 
into conservative and dissipative pieces, and review some 
properties of this decomposition derived by Mino [67| . 

We start by defining some notation. Suppose that we 
have a particle at a point V with four velocity m'', and 
that we are given a linearized metric perturbation h^^ 
which is a solution (not necessarily the retarded solution) 
of the linearized Einstein equation equation for which the 
source is a delta function on the geodesic determined by 
V and u^. The self-acceleration of the particle is then 
some functional of 7^, u^, h^n and of the spacetime metric 
5^1^, which we write as 

M- (2-69) 

Note that this functional does not depend on a choice of 
time orientation for the manifold, and also it is invari- 
ant under —u^. The retarded self-acceleration is 
defined as 

<ct [V. g^,] - a'' [V, g^,, h^;^] , (2.70) 

where is the retarded solution to the linearized Ein- 
stein equation obtained using the time orientation that 
is determined by demanding that be future directed. 
This is the physical self-acceleration which is denoted by 
a'' throughout the rest of this paper. Similarly, the ad- 
vanced self-acceleration is 

«adv I'P, "^ 5m^] - «^ [r, 5^., /if/] , (2.71) 

where is the advanced solution. It follows from these 
definitions that 

af^t [V, -u^ 5m-] = <dv "^ 5m-] ■ (2.72) 



We define the conservative and dissipative self- 
accelerations to be 

<ons - I «ct + <dv) > (2-73) 

and 

<^Ls = I i<t (2.74) 

The physical self-acceleration can then be decomposed as 

a^^a^,,^a^,,, + a'^,^,. (2.75) 

A similar decomposition applies to the forcing functions 
(EH: 

5i^' = 5il„s + 5i1iss: (2-76a) 
g['^ = Gil,. + Gil., (2.76b) 

for s = 1, 2. 

Next, we note that if ^ is any diffeomorphism from the 
spacetime to itself, then the self acceleration satisfies the 
covariance relation 

arct['4'{'P),ip*u'',ip*gf,^] ^ tp*a'^^^[V,u'',gf,^]. (2.77) 

Taking the point V to be {to,ro,9o,(j>Q) in Boyer- 
Lindquist coordinates, and choosing ip to he t —i- 2<o ~ t, 
(j) — + 200 ~ 4>i then ^ is an isometry, ^*g^u — 9iiv It 
follows that 

arot(-'"t,"r,Me, -u</>) = -ei^ajf^t (ut , , , M^) , (2.78) 
where 

£, = (1,-1,-1,1) (2.79) 

and there is no summation over v on the right hand side. 
Combining this with the identity (|2.72p gives 

al^^{ut,Ur, Me, Ucj,) = -e,aj:'j,t(ut, ~Ur, -ug, u^). (2.80) 

Now, under the transformation pr ~Pr, Pe ^ ^Pe 
with other quantities fixed, the action variables and the 
quantities Pa arc invariant, the angle variables q^. and qg 
transform as qr 2tt — qr, qg — » 2tt — qg, while qt — t and 
q,f, — (j) flip sign. This can be seen from the definitions 
([OT]) and (|2.35bp . Explicitly we have 

qtix'' , esps) - t = ~[qt{x'^,ps) ~t], (2.81a) 

q^{x'^ ,esps) - (f) = -[q4>ix'',ps) - (2.81b) 

qA{x'^,esPs) = 2tt ~ qA{x'^ ,ps), (2.81c) 

P,{x\esP5) = P^ix\ps), (2.81d) 

where we use the values (12.791) of £a, the functions qa are 
defined before Eq. I|2.57p . and qA = {qr,qe)- If we now 
differentiate with respect to pa holding x" fixed and use 
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the definitions ((237|) . (|2.53bp and (j2.59bp of the functions 
and we obtain 

/^(x^e^u^) = -e,mx^,u^), (2.82a) 
Fi^{x^,e^u^) = e^i^'^(a;'',M^). (2.82b) 



A similar computation gives 



We now compute the conservative and dissipative 
pieces of the forcing functions ^i^' and G\^\ using the 
definitions (PIM]) and Using the resuhs (p:5D|) 

and (|2.82p we obtain 



^!ret(^7'"7)' 



(2.84) 



J a advV"7 



K) = /aK)%adv("7) 

[-ei^/a(e7"7)] 

3irct(^7'^7)- 



i'"'i'rct(^7'"7) 



and using that the mapping x"^ x'^ , Ufj_ ^ 
sponds to Pj Pj, Qr ^ 271 — Qr, qg ^ 2tt 
(2.83) yields the identities 



finally 



and 



£Lns (QA .Pj) = gi'^ {qr ,q9,Pj)+ 9^'' (2^ - qr , 2^ 



si diss (9-4, Pj) 



G 



(1) 



^?dissilA. Pj) 



.gi^) {qr ,qe,Pj)- 9^^> (2^ - 9, , 2^ - ge , ) /2 



-,(1) 



{qA,Pj) = G['\qr,qe,Pj)-G['\2TT-qr,2TT-qg,Pj) /2, 



/2, 



Gi'\qr, qe, Pj) + G[^\27t - qr, 27T - qe, Pj) /2. 



(2.85a) 
(2.85b) 

(2.86a) 
(2.86b) 



Here we have used the fact that the forcing functions 
are independent of qt and g^, as discussed in the last 
subsection. Similar equations apply with gi^' and G,-"'^'' 

( s) ( s) 

replaced by the higher order forcing terms ga and G] , 
s > 2. 



variables 



It follows from the identity (|2.86ap that, for the action- 
variable forcing functions g[^\ the average over the 2- 
torus parameterized by qr and qe of the conservative piece 
vanishes. For generic orbits (for which ujr and LUg are in- 
commensurate) , the torus-average is equivalent to a time 
average, and so it follows that the time average vanishes, 
a result first derived by Mino [67]. Similarly from Eqs. 
(|2.85p it follows that the torus-average of the dissipative 
pieces of gi^^ vanish. 



q{t) = {qi{t),q2it),...,qN{t)), (3.1a) 
J(i) = (Ji(t), J2 (<),..., JA/(i)), (3.1b) 



and is defined by the equations 
dqa 



dt 
dJx 
dt 



cja(J, t) + £ga(q, J, t, e), 1 < a < A^, (3.2a) 



eG'A(q, J,?,e), 1 < A < M. 



(3.2b) 



Here the variable t is the "slow time" variable defined by 

i^et. (3.3) 



III. A GENERAL WEAKLY PERTURBED 
DYNAMICAL SYSTEM 



We assume that the functions ga and G\ can be expanded 
as 



In the remainder of this paper we will study in detail 
the behavior of a one-parameter family of dynamical sys- 
tems parameterized by a dimensionless parameter e. We 
shall be interested in the limiting behavior of the sys- 
tems as e ^ 0. The system contains N + M dynamical 



5a(q,J,f,e) = £5i'^)(q,J,tV-i 



gi'\ci,3,i)+g^^\ci,J,i)e + 0{e^), 

(3.4) 
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and 



= G['^ (q, J, t) + G\'' (q, J, t)s + 0{e'). 

(3.5) 

These series are assumed to be asymptotic series in e 
as e ^ that are uniform in We assume that the 
functions LOa, ga^ and G^^^ are smooth functions of their 
arguments, and that the frequencies cua are nowhere van- 
ishing. FinaUy the functions ga and Gx are assumed to 
be periodic in each variable qa with period 27r: 

5„(q + 2^k,J,i) = ga(q,J,t), 1 < a < A^, (3.6a) 
GA(q + 27rk,J,f) = GA(q,J,<), 1 < A < Af, (3.6b) 

where k = (fci, . . . , k^) is an arbitrary A^-tuple of inte- 
gers. 



(2), 



The equations (|2.47p derived in the previous sec- 
tion describing the inspiral of a point particle into 
a Kerr black hole are a special case of the dy- 
namical system (13.21) . This can be seen using the 
identifications t — q — (qt , Qr , qe , q4>) , J — 

iP2,P3,p4,M,,M2), Gi'^ = iG'i\G'i\Gi'\0,0) and 



(2) 



^ - {G^^\G'i\G''^\Gi,G2). The forcing functions 



G 

ga^ and G^"^ are periodic functions of qa since they de- 
pend only on the variables qA = {qr,Qe) which are angle 
variables; they do not depend on the variable qt which is 
not an angle variable. Note that the system (|3.2p allows 

the forcing functions yi*' , gI^^'* and frequencies uja to de- 
pend in an arbitrary way on the slow time t, whereas 
no such dependence is seen in the Kerr inspiral system 
(|2.47p . The system studied here is thus slightly more 
general than is required for our specific application. We 
include the dependence on t for greater generality and 
because it does not require any additional complexity in 
the analysis. 

Another special case of the system (|3.2I) is when N — 
M and when there exists a function H(3,i) such that 



Wa(J,i) 



dH{3, i) 
dJa 



(3.7) 



ior 1 < a < N. In this case the system p.2p represents 
a Hamiltonian system with slowly varying Hamiltonian 
H{3,t), with action angle variables {qa, Ja), and subject 



In other words, there exists T > such that for every q, J, every 
integer A'^, and every 5 > there exists €i = ei(q,J,Af, 5) such 



that 



JV 



(q, J, i, e)-Y, ffi"^ (q. J, 



<5e' 



to arbitrary weak perturbing forces that vary slowly with 
time. The perturbed system is not necessarily Hamilto- 
nian. 

Because of the periodicity conditions (|3.6p . we can 
without loss of generality interpret the variables qa to 
be coordinates on the A^-torus T^, and take the equa- 
tions (|3.2p to be defined on the product of this N-torus 
with an open set. This interpretation will useful below. 

In the next several sections we will study in detail the 
behavior of solutions of the system p.2p in the limit e 
using a two timescale expansion. We follow cl osely 
the exposition in the book by Kevorkian and Cole [7J], 
except that we generalize their analysis and also correct 
some errors (see Appendix [Bjl . For clarity we treat first, 
in Sec. IIVI the simple case of a single degree of freedom, 
N = M — 1. Section |V] treats the case of general A^ and 
AI, but with the restriction that the forcing functions 
ga and G\ contain no resonant pieces (this is defined in 
Sec. lVCp . The general case with resonances is treated in 
the forthcoming papers [ll, [zl]. Finally in Sec. IVll we 
present a numerical integration of a particular example 
of a dynamical system, in order to illustrate and validate 
the general theory of Sees. IIVI and fVl 



IV. SYSTEMS WITH A SINGLE DEGREE OF 
FREEDOM 

A. Overview 

For systems with a single degree of freedom the general 
equations of motion ()3.2p discussed in Sec. IIIII reduce to 

q{t) = uj{J,i) +eg{q,J,i,e), (4.1a) 
j{t) = eG{q,J,i,s), (4.1b) 

for some functions G and g, where t = st is the slow time 
variable. The asymptotic expansions p.4p and (|3.5p of 
the forcing functions reduce to 



l{q,J,i,e) = ^g(^)(g,J,tV-i 



s=l 



g<^'\q,J,i)+g(^\q,J,i)e + 0{e^), 



(4.2) 



and 



for all t with < t < T and for all e with < e < ei. 



G(g,J,f,£) = ^G(^)(g,J,tV-i 

= G(i)(g, J, i) + G(2) {q, J, i)e + Oie^). 

(4.3) 

Also the periodicity conditions (|3.6p reduce to 

g{q + 2n,J,i) = g{q,J,i), (4.4a) 
G{q + 271, J,i) = G{q,J,i). (4.4b) 
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In this section we apply two-timescale expansions to 
study classes of solutions of Eqs. I|4.ip in the limit e ^ 0. 
We start in Sec. IIVBI by defining our conventions and 
notations for Fourier decompositions of the perturbing 
forces. The heart of the method is the ansatz we make 
for the form of the solutions, which is given in Sec. IIV C] 
Sec. IIVDI summarizes the results we obtain at each or- 
der in the expansion, and the derivations are given in Sec. 
IIV El Although the results of this section are not directly 
applicable to the Kerr inspiral problem, the analysis of 
this section gives an introduction to the method of anal- 
ysis, and is considerably simpler than the multivariable 
case treated in Sec. FVl below. 

B. Fourier expansions of the perturbing forces 

The periodicity conditions (|4.4p apply at each order in 
the expansion in powers of e: 

g^'\q + 2^,J,t) = g^'\q,J,t), (4.5a) 
G'^'\q + 2Ti,J,i) = G^'Xq^J.i). (4.5b) 

It follows that these functions can be expanded as Fourier 
series: 

oo 

g(^)(g,J,i) = d^k^J^ty"^ (4-6a) 

k— — oo 
oc 

G^'\q,J,i) - Gi'\j,iy'''^ (4-6b) 

k— — oc 

where 

9k\j,t) = ^ r dqe-^'">g^^\q,J,i), (4.7a) 

Gk\j^t) = ^ I ^ dqe-^^^G^^\q,J,t). (4.7b) 

For any periodic function / = /(?), we introduce the 
notation 

{f) = ^ r f{q)dq (4.8) 
for the average part of /, and 

f{q) - f{q) - (/> (4.9) 

for the remaining part of /. It follows from these defini- 
tions that 

(5^^) (g, J, t)) = g(^) (J, i) , (G(^) {q, J, i)) = Gi'^ (J, i) , 

(4.10) 

and that 

9^'\q,J,i) = Ydk^M^''"' (4.11a) 
G(^)(g,J,?) = ^4^V,i>''^- (4.11b) 



We also have the identities 

iU - (/> = (4.12a) 
if 9) = {f9) + {.f){9) (4.12b) 

for any periodic functions f{q), g{q)- 

For any periodic function /, we also define a particular 
anti-derivative Jf of / by 

(X/)(g)^^^e^'=^ (4.13) 

where ft = J dqe~'^^'' f [q) / {2tt) are the Fourier coeffi- 
cients of /. This operator satisfies the identities 

(!/),« = /, (4.14a) 
{{If)g) = -{f{Ig)), (4.14b) 
(/(J/)) = 0. (4.14c) 

C. Two timescale ansatz for the solution 

We now discuss the ansatz we use for the form of the 
solutions of the equations of motion. This ansatz will be 
justified a posteriori order by order in e. The method 
used here is sometimes called the "method of strained 
coordinates" [t^ . 

We assume that q and J have asymptotic expansions 
in e as functions of two different variables, the slow time 
parameter t — et, and a phase variable ^' (also called 
a "fast-time parameter"), the dependence on which is 
periodic with period 2tt. Thus we assume 

oo 

q{t,e) = ^e^g(^)(vl/,t) 

s=0 

= q^"\^,i)+eq'^^H-i',i) + 0{e^), (4.15a) 

OO 

J{t,e) = 5]£V(^)(vI/,t) 

s=0 

= J(")(*,<) +ej(i)(*,i) + 0(e2).(4.15b) 

These asymptotic expansions are assumed to be uniform 
in i. The expansion coefficients J^-^^ are each periodic in 
the phase variable 5* with period 2tt: 

j('')(* + 27r,i) (*,?). (4.16) 

The phase variable 5* is chosen so that angle variable q 
increases by 2n when ^> increases by 27r; this implies that 
the expansion coefficients g'^^ satisfy 

g(°)(*-)-27r,i) = g(°)(*,i) + 27r, (4.17a) 
g(")(*-f 27r,i) = q'^'H^^f-.i), s > 1. (4.17b) 

The angular velocity il — d'i'/dt associated with the 
phase ^ is assumed to depend only on the slow time 
variable t (so it can vary slowly with time), and on e. 



20 



We assume that it has an asymptotic expansion in e as 
e — > which is uniform in t: 



'dt 



s=0 



(4.18) 
(4.19) 



Equation (|4.19p serves to define the phase variable ^ 
in terms the angular velocity variables Q^'^^i), s — 
0, 1, 2 . . ., up to constants of integration. One constant 
of integration arises at each order in e. Without loss 
of generality we choose these constants of integration so 
that 



g(^)(0,i) = 



(4.20) 



for all s, t. Note that this does not restrict the final 
solutions q{t,e) and J{t,e), as we show explicitly be- 
low, because there are additional constants of integra- 
tion that arise when solving for the functions g(*^(^,t) 
and J(^)(*,f). 

Roughly speaking, the meaning of these assumptions 
is the following. The solution of the equations of motion 
consists of a mapping from {t,s) to {q,J). That map- 
ping contains dynamics on two different timescales, the 
dynamical timescale ~ 1 and the slow timescale ~ 1/e. 
The mapping can be uniquely written the composition of 
two mappings 



(i,e) 



(q, J), 



(4.21) 



such that the first mapping contains all the fast dynam- 
ics, and is characterized by the slowly evolving frequency 
n{t,e), and the second mapping contains dynamics only 
on the slow timescale. 



D. Results of the two-timescale analysis 



By substituting the ansatz (j4.15bp - (|4.20p into the 
equations of motion (|4.1[) we find that all of the assump- 
tions made in the ansatz can be satisfied, and that all of 
the expansion coefficients are uniquely determined, or- 
der by order in e. This derivation is given in Sec. IIV El 
below. Here we list the results obtained for the various 
expansion coefficients up to the leading and sub-leading 
orders. 



1. Terminology for various orders of the approximation 

We can combine the definitions just summarized to 
obtain an explicit expansion for the quantity of most in- 
terest, the angle variable g as a function of time. From 
the periodicity condition (|4.17ap it follows that the func- 
tion (vj/, t) can be written as -I- g'"' i) where g^"' 
is a periodic function of [We shall see that q^^^^ in fact 



vanishes, cf. Eq. (|4.27p below.] From the definitions (|3.3p 
and (|4.19p . we can write the phase variable 5* as 

* = (t) + 7/^(1) (<) + e^A^') {i) + 0{e^), (4.22) 

e 

where the functions tp^'^^ (i) are defined by 

Inserting this into the expansion (|4.15aP of q and using 
the above expression for g'"' gives 



it,e) = i^(°)(i)+ [V'^'Hi) + 9^°^(*,?) 
e L 

+ g(i)(5-,t)j + 0(e2). (4.24) 



1 

e 

+e 



We will caU the leading order, 0(l/e) term in Eq. (|4.24p 
the adiabatic approximation, the sub- leading 0(1) term 
the post- 1- adiabatic term, the next 0{e) term the post-2- 
adiabatic term, etc. This choice of terminology is moti- 
vated by the terminology used in post-Newtonian theory. 

It is important to note that the expansion in powers 
of £ in Eq. (14. 24^ is not a straightforward power series 
expansion at fixed t, since the variable ^' depends on e. 
[The precise definition of the expansion of the solution 
which we are using is given by Eqs. (|4.15ap - (|4.20p .] 
Nevertheless, the expansion (|4.24p as written correctly 
captures the e dependence of the secular pieces of the 
solution, since the functions q^°^ and q^'^'' are periodic 
functions of and so have no secular pieces. 



2. Adiabatic Order 

First, the zeroth order action variable is given by 

j(o)(^,f) = J7(0)(f), (4.25) 
where J'^^'' satisfies the differential equation 
dj(°^{i) 



dt 



= G(^)[^W(f),t]. 



(4.26) 



Here the right hand side denotes the average over q of 
the forcing term G'^^^q, cf. Eqs. dM]) above. 

The zeroth order angle variable is given by 



g(")(vl/,f) = *, 



(4.27) 



and the angular velocity Q that defines the phase variable 
^' is given to zeroth order by 



n(°\i) = uj[j^°Hi),t\. 



(4.28) 



Note that this approximation is equivalent to the follow- 
ing simple prescription: (i) Truncate the equations of 
motion (j4.ip to the leading order in e: 

q{t) = uj{J,i) + eg^-^\q,J,t), (4.29a) 
i(t) = eG^^\q,J,i)- (4.29b) 
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(ii) Omit the driving term g^^^ in the equation for the 
angle variable; and (iii) Replace the driving term G*-^^ in 
the equation for the action variable with its average over 

q- 



3. Post-l-adiabatic Order 



Next, the first order action variable is given by 

,-,^2SVHM + j(..,£), (430, 

where the symbol X on the right hand side denotes the in- 
tegration operator l|4.13p with respect to In Eq. (|4.30p 
the quantity (i) satisfies the differential equation 



dt 



dJ 



(G(i).9 



17(0) (t) f^(0)(f) 



G 



(2) 



(4.31) 



Here it is understood that the quantities on the right 
hand side are evaluated at q = q^'^^ = ^ and J — J^'^\i). 
The sub-leading correction to the phase variable 5" is 
given by 

= ^£j[J^'\hi]J^'\t)+9'^^\j^°\i)A (4.32) 
Finally, the sub-leading term in the angle variable is 
g(i)(*,t) = q^^\^,i) + Q(i)(t), (4.33) 

where 

1 duj, 



and 



QW(t) = -g«(0,t). 



4-. Discussion 



(4.34) 



(4.35) 



One of the key results of the general analysis of this 
section is the identification of which pieces of the external 
forces are required to compute the adiabatic and post-l- 
adiabatic solutions. From Eqs. (g^ll), and ([i?^ . 
the adiabatic solution depends only on the averaged piece 
Gq"^' (J, i) — (G'-^) (g, J, t)) of the leading order external 
force G^^). This quantity is purely dissipative, as can be 
seen in the Kerr inspiral context from Eqs. (|2.86p and 
(|2.85p . More generally, if the perturbing forces g and 



G arise from a perturbation e/S.H = E^Aff*^") to the 
Hamiltonian, then the forcing function G'^-* is 



G(^)((Z,J,i) 



dq 



and it follows that the average over q of G^*-* vanishes. 

At the next order, the post-l-adiabatic term 
depends on the averaged piece G^\j, i) — (G^-^^ (g, J, i)) 
of the sub- leading force G^^^ , again purely dissipative, as 
well as the remaining conservative and dissipative pieces 
of the leading order forces G'^^^{q,J,i) and g^^\q, J,i). 
This can be seen from Eqs. (|4.3ip and (|4.32p . These re- 
sults have been previously discussed briefly in the EMRI 
context in Refs. [13, [11]. For circular, equatorial orbits, 
the fact that there is a post-l-adiabatic order contribu- 
tion from the second order self-force was first argued by 
Burko [H. 



5. Initial conditions and the generality of our ansatz 

We will show in the next subsection that our ansatz 
(|4.15ap - (|4.20p is compatible with the one parameter 
family of differential equations (|4.ip . However, it does 
not necessarily follow that our ansatz is compatible with 
the most general one parameter family [q{t, e), J{t, e)] of 
solutions, because of the possibility of choosing arbitrary, 
e-dependent initial conditions g(0,e) and J{0,e) at the 
initial time t = 0.^^ In general, the e dependence of 
the solutions arises from both the e dependence of the 
initial conditions and the e dependence of the differential 
equations. It is possible to choose initial conditions which 
are incompatible with our ansatz. 

To see this explicitly, we evaluate the expansions (|4.24p 
and ((430)) at i i = 0. This gives 



g(0,e) 
J(0,e) 



.-i^(o)(o) + ^(i)(0) + 0(e). 



(4.36a) 



xG(i' [e- (0) + V'(i) (0), (0), 0] 



(4.36b) 



Recalling that parameters ip'-^HO), ip'-^^O), j'-°'>{0) and 
J^^^ (0) are assumed to be independent of e, we see that 
the conditions (|4.36p strongly constrain the allowed e de- 
pendence of the initial conditions. We note, however, 
that the choice of constant (e independent) initial condi- 
tions 



g(0,e)=qo, J(0,ir) = Jo 



(4.37) 



More generally we could consider specifying initial conditions at 
some time t = tg. In that case we would modify the definition of 
the rescaled time coordinate to i = e{t — to). 
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can be accommodated, which is sufBcient for most appH- 
cations of the formahsm. To achieve this one chooses 



4'^°\0)^0, V^i)(0) = go, J^'\0) = Jo, (4.38) 



and 



XG(i)[go, Jo,0] 
^[^0,0] 



(4.39) 



Derivation 



In this subsection we give the derivation of the results 
(|4.25p - (|4.35p summarized above. At each order s we 
introduce the notation J7^''^(t) for the average part of 
J('')(*,i): 



J^'\i) = (J(")(*,t)) = 7^ / J(")(*,f)d*. (4.40) 
27r Jo 

We denote by J^*' the remaining part of J^^\ as in Eq. 
(14. 9p . This gives the decomposition 

J^'H*,^) = J^'Hi) + j'-"H'^,i) (4.41) 

for aU s > 0. Similarly for the angle variable we have the 
decomposition 



(4.42) 



for all s > 1. [We do not use this notation for the s ~ 
case for the angle variable, since is not a periodic 
function of ^f, by Eq. (|4.17ap ]. 

Using the expansions (|4.15ap and (|4.15bp of q and J 
together with the expansion (|4.19p of d'^/dt, we obtain 



dq 
'dt 



(0) 



+Oie^). 



(4.43) 



Here we use commas to denote partial derivatives. We 
now insert this expansion together with a similar expan- 
sion for dJ /dt into the equations of motion (|4.ip and use 
the expansions (|4.2p and (|4.3p of the external forces g 
and G. Equating coefficients^^ of powers of e then gives 
at zeroth order 



(4.44a) 
(4.44b) 



As is well known, this procedure is valid for asymptotic series as 
well as normal power series. 



at first order 

~ c.,^j(i) = -f^^i)?!^^ - qf + .g(^\(4.45a) 
^(1) j(0) j^(o) ^ _ j(0) ^ gd) ^ (4 45b) 

and at second order 

+5(^)-J7(^)g(^)-f7(i)gW 
-gj\ (4.46a) 

-//^+G(2). (4.46b) 

Here it is understood that all functions of g and J are 
evaluated at g'"' and j'"' . 

1. Zeroth order analysis 



The zeroth order equations (j4.44p can be written more 
explicitly as 

f](0)(t)g;:^)(vl/,t) = c.[j(")(vl/,i),t], (4.47a) 
n^''\t)J%\-^,i) = 0. (4.47b) 

The second of these equations implies that J^^-' is inde- 
pendent of so we obtain J^^\^ ,i) = J^°\i). The 

first equation then implies that g^^^ is independent of 
and integrating with respect to ^' gives 



* + s(°)(f), 



(4.48) 



where Q*-"-* is some function of i. The periodicity con- 
dition (|4.17ap now implies that the coefficient of ^' in 
Eq. (|4.48p must be unity, which gives the formula (|4.28p 
for the angular velocity ri(°)(t). Finally, the assumption 
()4.20p forces Q'-*'-'(t) to vanish, and we recover the for- 
mula 14:27)1 for g(°)(*,i). 



2. First order analysis 



The first order equation (j4.45bp can be written more 
explicitly as 

f}(°)(i)j«(vl',t) = ^jfii) 

-KG(i)[*,^(")(f),t], (4.49) 

where we have simplified using the zeroth order solutions 
(|4.25p and (|4.27p . We now take the average with respect 
to of this equation. The left hand side vanishes since 
it is a total derivative, and we obtain using the definition 
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(|4J|) the differential equation (14261) for Next, 
we subtract from Eq. (|4.49p its averaged part, and use 
the decomposition ()4.4ip of J*^^) . This gives 

r!(")(<)/^)(*,<) = G(i)[*,:7(")(t),i]. (4.50) 

We solve this equation using the Fourier decomposition 
(|4llbl) of G(i) to obtain 



fe#0 



(4.51) 



This yields the first term in the result (|4.30p for j'-'^^ when 
we use the notation (|4.13p . 

Next, we simplify the first order equation (|4.45ap using 
the zeroth order solutions (j4.25p and (|4.27[) . to obtain 

{i)q%^ i) - (t), i]J<^'^ i] 

= -r!«(i) + g(^^[^,J^"Hi),i]. (4.52) 

Averaging with respect to 5" and using the decomposi- 
tions ()4.4H) and (I4.42p of J^^^ and g^^^ now gives the 
formula (1432]) for ri(^)(t). Note however that the func- 
tion J'^^^ (t) in that formula has not yet been determined; 
it will be necessary to go to one higher order to compute 
this function. 

Finally, we subtract from Eq. (|4.52p its average over 
^ using the decompositions (j4.4ip and (|4.42p and then 
integrate with respect to '5 using the notation (|4.13p . 
This gives 



c.,j[j(°Ht),t]Jj(i)[vI/,t] 



f7(")(t) 



(4.53) 



Using the result for J^^^ given by the first term in Eq. 
(|4.30p now yields the formula (|4.34p for g^^) (*,?), and the 
result (|4.33p for g*^^) then follows from the decomposition 
(|4.42p together with the initial condition (|4.20p . 



3. Second order analysis 



We simplify the second order equation (j4.46bp using 
the zeroth order solutions (|4.25p and (|4.27p , average over 
^, and simplify using the decompositions (|4.4ip and 
(|02)) and the identities (|4?T2l) . The result is 

+ (g«(*,i)G«[vI/,^(o)(i),t]; 



+ [J^'>i^, i) G[]>[<f, J^''Ki)A) ■ (4.54) 

Using the expressions (|4.34p and (|4.30p for q^^^ and j'-'^' 
and simplifying using the identities (|4.14p now gives the 
differential equation (|4.3ip for j/'^^. 



^. Extension to arbitrary order 

In this subsection we prove by induction that solutions 
are uniquely determined at each order in e. Our inductive 
hypothesis is that, given the equations up to order s, we 
can compute all of the expansion coefficients g*^"^(4',t), 
j(")(*,i) and r2(")(t) for < m < s, except for the av- 
eraged piece J^'\i) of J(")(*,f), and except for Vl'^'^t), 
which is assumed to be determined by j'^'^'^ (t) . From the 
preceding subsections this hypothesis is true for s = 
and for s = 1. We shall assume it is true at order s — 1 
and prove it is true at order s. 

The equations of motion at order s, when simplified 
using the zeroth zeroth order solutions (j4.25p and (|4.27p . 
can be written as 



+S (4.55a) 



(4.55b) 



Here S = 5(^',t) and T — T{'i>,i) are expressions in- 
volving the forces G'"-* and g*^"' for < u < s evaluated 
at g = = and J = = and involving the 

coefficients and for < u < s - 2 which 

by the inductive hypothesis are known. Therefore we can 
treat S and T as known functions. 

Averaging Eq. (|4.55bp over yields the differential 
equation 

+ {g[]U^'-''>). (4.56) 

By the inductive hypothesis all the terms on the right 
hand side are known, so we can solve this differential 
equation to determine J^^^~^\ 

Next, averaging Eq. (|4.55ap yields 

{S). (4.57) 



(s-l) 



Since J''-*^^'' has already been determined, the right hand 
side of this equation is known and therefore the equation 
can be used to solve for Q^'^'> once J'^'^^ is specified, in 
accord with the inductive hypothesis. Next, Eq. (j4.55bp 
with the average part subtracted can be used to solve 
for J^'^^ and once J^^^ is known Eq. (j4.55ap with the 
average part subtracted can be used to solve for g^*^. 
Finally, the averaged piece can be 
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computed from g*^^) using the initial condition (|4.20p and 
the decomposition (I4.42p . Thus the inductive hypothesis 
is true at order s if it is true at order s — 1 . 



B. Fourier expansions of perturbing forces 



V. SYSTEMS WITH SEVERAL DEGREES OF 
FREEDOM SUBJECT TO NON-RESONANT 
FORCING 



Overview 



In this section we generahze the analysis of the preced- 
ing section to the general system of equations (|3.2p with 
several degrees of freedom. For convenience we reproduce 
those equations here: 



dqg 
dt 

dJx 
dt 



Ljg (J, i) + eg'i^ (q, J, f) + e2g(^2) 

+0{e^), l<a<N, (5.1a) 

eGi'^(q,J,f)+e2G'f (q,j,f) 



+0{e^), 1 < A < M. 



(5.1b) 



For the remainder of this paper, unless otherwise spec- 
ified, indices a,/?, 7,(5, e, . . . from the start of the Greek 
alphabet will run over 1 . . . , and indices A, /i, v, p,a, . . . 
from the second half of the alphabet will run over 1 . . . M. 

The generalization from one to several variables is 
straightforward except for the treatment of resonances 
[T^I- The key new feature in the TV variable case is that 
the asymptotic expansions now have additional terms 
proportional to ^/e, e"^/^, ... as well as the integer pow- 
ers of e. The coefficients of these half-integer powers of 
£ obey source- free differential equations, except at reso- 
nances. Therefore, in the absence of resonances, all of 
these coefficients can be set to zero without loss of gen- 
erality. In this paper we develop the general theory with 
both types of terms present, but we specialize to the case 
where no resonances occur. Subsequent papers [zl, [t^ 
will extend the treatment to include resonances, and de- 
rive the form of the source terms for the half-integer 
power coefficients. 

We start in Sec. IV Bl by defining our conventions and 
notations for Fourier decompositions of the perturbing 
forces. In Sec. IV Cl we discuss the assumptions we make 
that prevent the occurrence of resonances in the solu- 
tions. The heart of the method is the ansatz we make 
for the form of the solutions, which is given in Sec. IV Dl 
Section IV El summarizes the results we obtain at each 
order in the expansion, and the derivations are given in 
Sec. IV Fl The implications of the results are discussed in 
detail in Sec. IVIII below. 



The periodicity condition (|3.6[) applies at each order in 
the expansion in powers of e, so we obtain 

gi'\q + 27Tk,J,i) = gi^)(q,J,f), (5.2a) 
G^^\q + 27Tk,3,i) = G^^\<i,J,i), (5.2b) 

for s > 1, 1 < a < iV, and 1 < A < M. Here 
k — (fci, . . . , /cjv) can be an arbitrary A^-tuple of inte- 
gers. It follows from Eqs. (|5.2p that these functions can 
be expanded as multiple Fourier series: 

5i^)(q,J,f) = E-9il(J'^~)^'''^ (5-3a) 



Gi'^)(q,J,f) = ^Gil(J,f)e*-' 



(5.3b) 



where 



(J,<) = y'd^ge-*'^-'i5i^)(q,J,f), (5.4a) 

G^xli^J) = y'd^ge-'^ 'iGi'^Hq,J,0.(5.4b) 

Here we adopt the usual notations 

oo oo 



(5.5) 



k ki— — oo kN— — oo 



d''q= I dqi... dqN. (5.6) 
Jo 



and 



N 



k • q = ^ kaqa- 



(5.7) 



For any multiply periodic function / — /(q), we intro- 
duce the notation 



(/> = 



1 



d^qfiq) 



(27r)^ 

for the average part of /, and 

/(q) - /(q) - (/) 



(5.8) 



(5.9) 



for the remaining part of /. It follows from these defini- 
tions that 

(gi^)(q,J,f)) = .9i'^^(J,f), (Gi^)(q,J,f)) =Git(J,f), 

(5.10) 

5i^)(q,J,f) = E5il(J>^>^''-^ (5-lla) 
G(^)(q,J,f) = Y.G^^l{3,i)e^-^. (5.11b) 



and that 
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We also have the identities 
,dqa 

ifg) = (/3> + (/>(5> 



(5.12a) 
(5.12b) 



for any multiply periodic functions /(q), 9i^)■ 

For any multiply periodic function / and for any vector 
V = (wi, . . . , wjv), we also define the quantity J^f by 



(^v/)(q) - E TIT 



(5.13) 



where fk — J d^qe ^^''^f{q)/{2TT)^ are the Fourier co- 
efficients of /. The operator 1^ satisfies the identities 



2'v(v-V/) = /, (5.14a) 
((Xv/)g) = -(/(Xvff)), (5.14b) 
(/(Xv/)) - 0. (5.14c) 

C. The no- resonance assumption 

For each set of action variables J and for each time t, 
we will say that an N-tuple of integers k 7^ is a resonant 
N-tuple if 



k-a;(J,t) =0. 



(5.15) 



where uj = (lji, . . . ,lun) are the frequencies that appear 
on the right hand side of the equation of motion (j3.2ap . 
This condition governs the occurrence of resonances in 
our perturbation expansion, as is well known in context of 
perturbations of multiply periodic Hamiltonian systems 
|95| . We will assume that for a given k, the set of values 
of t at which the quantity 



(5.16) 



vanishes (i.e. the resonant values) consists of isolated 
points. Here is the leading order solution for 

J given by Eq. (|5.29p below. This assumption excludes 
persistent resonances that last for a finite interval in t. 
Generically we expect this to be true because of the time 
dependence of J'^^^'it). 

Our no-resonance assumption is essentially that the 
Fourier components of the forcing terms vanish for res- 
onant N-tuples. More precisely, for each fixed k and for 
each time t-c for which CTk(ir) = 0, we assume that 



9aL 



G 



is) 
Ak 



= 0, 
= 0, 



(5.17a) 
(5.17b) 



for s > 1 and for all t in a neighborhood of t,-. Our no- 
resonance assumption will be relaxed in the forthcoming 
papers 



In our application to inspirals in Kerr black holes, the 
no-resonance condition will be automatically satisfied for 
two classes of orbits: circular and equatorial orbits. This 
is because for these orbits there is either no radial mo- 
tion, or no motion in 0, and so the two-dimensional torus 
{qr,qe) reduces to a one-dimensional circle. The reso- 
nance condition k^LOr + kgujQ = reduces to kj-LOr = 
for equatorial orbits, or kgug = for circular orbits, and 
these conditions can never be satisfied since the funda- 
mental frequencies LOr and loq are positive. 



D. Two timescale ansatz for the solution 

We now discuss the two-timescale ansatz we use for the 
form of the solutions of the equations of motion. This 
ansatz will be justified a posteriori order by order in ^/e. 
Our ansatz essentially consists of the assumption that 
the mapping from [t, s) to (q, J) can be written as an 
asymptotic expansion in ^/e, each term of which is the 
composition of two maps, the first from {t, s) to an ab- 
stract iV-torus with coordinates = (^'i, . . . , ^'at), and 
the second from {^,t,e) to (q, J). Here t = et is the 
slow time parameter. All the fast timescale dynamics is 
encapsulated in the first mapping. More precisely, we 
assume 



(*, i) + (*, i) + (*, i) 

+eV\(y^){^J) + 0{e^), (5.18a) 



Mt,s) = E^"^'-^i"^'^(*'*~) 



/I' 



jf ^ (*, i) + V~eJi'^'^ (*, i) + eji'^ (*, i) 
+e3/2jf/'^(*,f) + 0(e2). (5.18b) 



These asymptotic expansions are assumed to be uni- 

(s) 

form in t. The expansion coefficients , where s = 
0, 1/2, 1, . . ., are multiply periodic in the phase variables 
with period 27r in each variable: 



j(^)(* + 27rk,<) = j(^)(*,t). 



(5.19) 



Here k — (fci, . . . , kj^) is an arbitrary A^-tuple of integers. 
The mapping of the abstract iV-torus with coordinates 
^ into the torus in phase space parameterized by q is 
assumed to have a trivial wrapping, so that the angle 
variable qa increases by 27r when increases by 2tt; 
this implies that the expansion coefficients q*^"^ satisfy 

gi°)(* + 2^k,i) - ql°\^,i) + 2Trk^, (5.20a) 
qi'H^ + 2TTk,i) = q^:\^,i), s > 1/2, (5.20b) 

for arbitrary k. The variables ^'i, . . . , ^'jv are sometimes 
called "fast-time parameters" . 
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The angular velocity 



where the functions 'il}a\i) are defined by 



(5.21) 



associated with the phase ^'^ is assumed to depend only 
on the slow time variable t (so it can vary slowly with 
time), and on e. We assume that it has an asymptotic 
expansion in ^/e as e ^ which is uniform in t: 



(5.22) 



71=0 



+e3/2fi(3/2)(i) + 0(£2). (5.23) 



Equations (|5.2ip and (|5.23p serve to define the phase vari- 
able ^'q, in terms the angular velocity variables Qa\i), 
s — 0, 1/2, 1 . . ., up to constants of integration. One con- 
stant of integration arises at each order in y^, for each 
a. Without loss of generality we choose these constants 
of integration so that 



(5.24) 



for all a, s and i. Note that this does not restrict the 
final solutions qa{t, e) and J\{t, e), as we show explicitly 
below, because there are additional constants of integra- 
tion that arise when solving for the functions gi'*'' ('S', t) 
(*,?). 



and 



(5.26) 



l^(°)(f) + i=V.i^/^)(?) 



Inserting this into the expansion (|5.18ap of qa gives 
qa{t,e) 

t V 

+V^[p^^/'Hi) + qi'^'H^,i) 

+8 Ui')(f) + g«(*,i)l + 0{e^/'). (5.27) 



We will call the leading order, 0(l/e) term in Eq. (|5.27p 
the adiabatic approximation, the sub-leading 0(1 /\/e) 
term the post- 1/2- adiabatic term, the next 0(1) term 
the post- 1- adiabatic term, etc. Below we will see that the 

functions qi"'' and qa^'^"' in fact vanish identically, and 
so the oscillatory, 'J'-dependent terms in the expansion 
(|5.27p arise only at post-2-adiabatic and higher orders. 

As before we note that the expansion in powers of e in 
Eq. (|5.27p is not a straightforward power series expansion 
at fixed i, since the variable ^' depends on e. [The precise 
definition of the expansion of the solution which we are 
using is given by Eqs. I|5.18ap - (|5.24p .] Nevertheless, 
the expansion (|5.27p as written correctly captures the e 
dependence of the secular pieces of the solution, since 
the functions q^°\ qa^"^^ and qa^ are multiply periodic 
functions of 'i' and so have no secular pieces. 



E. Results of the two-timescale analysis 

By substituting the ansatz (|5.18bp - (|5.24p into the 
equations of motion (|3.2p we find that all of the assump- 
tions made in the ansatz can be satisfied, and that all of 
the expansion coefficients are uniquely determined, or- 
der by order in ^/£. This derivation is given in Sec. IV Fl 
below. Here we list the results obtained for the various 
expansion coefficients up to the first three orders. 



1. Terminology for various orders of the approximation 

We can combine the definitions just summarized to 
obtain an explicit expansion for the quantity of most in- 
terest, the angle variables qa as a function of time. From 
the periodicity condition (|4.17ap it follows that the func- 
tion qa\^,i) can be written as ^'q, -I- qa\^,i) where 

qa^ is a multiply periodic function of '5'. From the defi- 
nitions (|3.3p and (|5.23p , we can write the phase variables 
^'q as 

*a = ^ V'i"^ (?) + ^^i^/^^ {i) + ^i^) (i) + Vi^i^/^) it) 

+eVi'Hi) + 0(£3/2), (5.25) 



2. Adiabatic Order 

The zeroth order action variables are given by 

jf (*,t) = :7f (t), (5.28) 

where J^°\i) = . . .,J^\i)) satisfies the set of 

coupled ordinary differential equations 




= G«[J^(°)(i),t]. (5.29) 



Here the right hand side denotes the average over q of 
the forcing term G^^^q, J^°'> cf. Eqs. (EH) above. 
The zeroth order angle variables are given by 

q^°\^,i) = ^a, (5.30) 

and the angular velocity Qa that defines the phase vari- 
able '^a is given to zeroth order by 

n^^Hi)=u;a[J^°\i),i]. (5.31) 

Note that this approximation is equivalent to the fol- 
lowing simple prescription: (i) Truncate the equations of 
motion (|5.ip to the 0{e); (ii) Omit the driving terms gi^' 
in the equations for the angle variables; and (iii) Replace 
the driving terms G^^^ in the equations for the action 
variables with their averages over q. 
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3. Post-l/2-adiahatic order 



Next, the 0{^) action variables are given by 

= (5.32) 

where J^^l''\t) = . . . , Ju'^Kt)) satisfies the 

set of coupled, source-free ordinary differential equations 



^ ' ^[J^'\i)M'/'\i) = 0. (5.33) 



.(1 



dt 



dJ, 



Equation (|5.33p will acquire a source term in Ref. [79| 
where we include the effects of resonances. The 0{y/e) 
angle variables are given by 



^/2)(*,t)=0, 



(5.34) 



and the angular velocity fla that defines the phase vari- 
able is given to 0{^/e) by 

^i'^'Hi) = ^[J^'\i)M'^"\i)- (5.35) 



Note that Eqs. (|5.33p and (|5.35p can be obtained sim- 
ply by linearizing Eqs. (|5.29[) and (|5.31[) about the ze- 
roth order solution. That is, J"^"^ + y/ej'^'^^^^ and 
rj'"' + VeS^*^^^' satisfy the zeroth order equations ([5?29l) 
and ((53T|) to 0{^/e). This means that setting J-(i/2) 
and rj(i/2) 

zero does not cause any loss of generality 
in the solutions (under the no-resonance assumption of 
this paper), as long as we allow initial conditions to have 
sufhciently general dependence on e. 



4- Post-l-adiabatic order 

The first order action variable is given by 

(*,tO =Xj^<„,^^-^G«[*, J^(")(t),i]+ J'«(i), (5.36) 

where the symbol I on the right hand side denotes the 
integration operator (|5.13p with respect to '4', gI^^^ is the 
non-constant piece of G^^"* as defined in Eq. (|5.9p . and 
is given by Eq. ((OT|) . In Eq. ((06l) the quantity 
JT"*-^^ {i) satisfies the differential equation 



di dJn 



(^(2) , ^ O ^AO jil/2) j{l/2) 



(5.37) 



Here it is understood that the quantities on the right 
hand side are evaluated at J = and q = q^*^^ = 

^. The last three terms on the right hand side of Eq. 
(|5.37p can be written more explicitly using the definition 
(|5.13p of T and the definition (jS.Sp of the averaging (...) 
as 



ik. 



^(i)*^(i) 



^ o(o) 



k G^^^*o^^^ 

'^""^Ak .Vak 



dJ„ 



(5.38) 



by 



The 0(e) correction to the angular velocity is given 



1 d'^LJa 

'2dJxdJ, 



(5.39) 

Finally, the sub-leading term in the angle variable is 



gW(vI/,f)=gli)(*,i) + QW(t), (5.40) 



where 



duJa 

'dTx 



[jr(")(t),t] 



+%(o)(,~)5i'^[*,.7(°) (5.41) 

QL'H?) = -9i'H0,f). (5.42) 
5. Discussion 



and 



One of the key results of the general analysis of this 
section is the identification of which pieces of the exter- 
nal forces are required to compute the adiabatic, post- 
1/2-adiabatic and post-l-adiabatic solutions. From Eqs. 
(|5.29p . (|5.3ip and (|5.27p . the adiabatic solution depends 

only on the averaged piece G^y^]^(i,i) = (G^^^ (q, J, f)) of 

the leading order external force G^^'' . Only the dissipative 

piece of the force G^^^ normally contributes to this aver- 
age. For our application to inspirals in Kerr, this follows 
from the identity (|2.86ap which shows that the average 

of the conservative piece of G^^"* vanishes. For a Hamil- 
tonian system with N — M, if the perturbing forces ga 
and Gp arise from a perturbation eAH = J^s e^Ai? 



to the Hamiltonian, then the forcing function G^ is 



aAi?("'(q, J,f) 
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and it follows that the average over q of G^*-* vanishes. 

At the next, post-l/2-adiabatic order, it follows from 
Eqs. and ([OS]) that the term ip^^^\i) depends 

again only on the averaged, dissipative piece G^^q of the 
leading order force. However, we shall see in the forth- 
coming paper [79| that when the effects of resonances 
are included, additional dependencies on the remaining 
(non-averaged) pieces of the first order self forces will 
arise. 

At the next, post-l-adiabatic order, the term il)a\i) 
in Eq. (j5.27p depends on the averaged piece G\Q{3,t) = 

{G^^\(l,3,i)) of the sub-leading force G^^\ again nor- 
mally purely dissipative, as well as the remaining conser- 
vative and dissipative pieces of the leading order forces 
gI^^-* (q, J, i) and ^ (q, 3,i). This can be seen from Eqs. 
(I5.37P and (|5.39p . These results have been previously 
discussed briefly in the EMRI context in Refs. [13, [68| . 
For circular, equatorial orbits, the fact that there is a 
post-l-adiabatic order contribution from the second or- 
der self- force was first argued by Burko (soj . 

Finally, we consider the choice of initial conditions for 
the approximate differential equations we have derived. 
The discussion and conclusions here parallel those in the 
single variable case, given in Sec. IIVD 51 above, and the 
results are summarized in Sec. IVII CI below 



where q^\'^,i) is multiply periodic in '5' with zero av- 
erage. 

Using the expansions (|5.18ap and (|5.18b[) of Qa and 
together with the expansion (|5.23p of d'i'a/dt, we obtain 



dqa 
dt 



(1/2)^(0) 



,(0)„(l/2) 



,(1/2)^(1/2) 



(0)„(1) 



JO) 



3/2 



f7(3/2)J0) +r!^l)a(l/2) , f^(l/2) (1) 



+ , „(l/2) 



"0 



^(3/2) (1/2) ^(1) (1) f^(l/2)^(3/2) 



+"/3 ya,f^ 



(5.48) 



We now insert this expansion together with a similar ex- 
pansion for dJx/dt into the equations of motion (|3.2p and 
use the expansions p.4p and (|3.5p of the external forces 
and G\. Equating coefhcients of powers^'^ of ^/e then 
gives at zeroth order 



"/3 ya,*3 — 

o(o) 7(0) - 

"/3 "^A,*3 ~ 



0, 



(5.49a) 
(5.49b) 



at order 0{\/e) 



F. Derivation 

We will denote by H-it) the set of resonant N-tuples 
k at time t, and by TZ'^(t) the remaining non-resonant 
nonzero N-tuples. The set of all N-tuples can therefore 
be written as the disjoint union 



{o}u7^(^)u7^"(^). 



(5.43) 



At each order s we introduce the notation J'jf'' (i) for 
the average part of J^^\'^,i): 



(5.44) 



(27r)^ Jo 



2n 



d-^l 



In 



d*Arj(^)(*,t). 



We denote by the remaining part of , as in Eq. 
(|5.9p . This gives the decomposition 

j'-^\^:t)^j^\t)^j't\^:t) (5.45) 

for all s > 0. Similarly for the angle variable we have the 
decomposition 

(*, t) = Qi'^) (t) + <jt ) (*, f) (5.46) 

for all s > 1/2. For the case s = we use the fact that 
qa\'*i^,i) — is a multiply periodic fimction of by 
Eq. (I5.20ap . to obtain the decomposition 

g^o) (*, i) = + Q(^o) ^ ^(0) ^~)^ (5^47) 



^(0) (1/2) 

o(0) r(i/2) 

"/3 •^A.l'a 



^(1/2) (0) .(1/2) 
0(l/2) t(0) 



at order 0{e) 

q(0)„(1) _ -0(l/2) (1/2) _ Q^^^flW 



(i)„(o) 



-,(0) 



(5.50a) 
(5.50b) 

^(1) 



+c.„,,,jl^' + ic.„,,,,,jl^/^'j(V2),(5.51a) 



0(0) t(1) _ _o(l/2) 7(1/2) 
"/3 "^A,*^ ~ "/? ""^ 



(1) 7(") 



13 '^A,*^j 



7(0) 
A,t 



-HG 



(1) 

A ' 



(5.51b) 



at order 0{e^^'^) 

o(0)„(3/2) _ 
"/3 ya,*^ ~ 



(-,(1/2) (1) o(l)«(l/2) o(3/2)„(0) 

+ i-„,,,,^,„ji^/^)j(V2)j(V2), (5.52a) 



0(0) 7(3/2) _ 

"/3 "^A,*^ — 



0(1/2) .(1) _ ^{1) .(1/2) 
" A 



/3 ''A.t/j "/3 '^A,*^ "/3 
7(1/2) ^(1) (1/2) ^(1) (1/2) 

(5.52b) 



This is justified since both sides are asymptotic expansions in 
powers of yTs at fixed i. 
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and at order O(e^) 



"/3 la,'!' 13 



o(o) 7(2) 



_(-,(l/2) (3/2) _ o(l)„(l) _ o(3/2)„(l/2) 
,1(1) t(1/2) .(1/2) (1) (1/2) .(1/2) 

+LU 7 I^^^ + -UJ , T T 7(1)7(1/2)7(1/2) 

+ —UJ T r T T 7(1/2) 7(1/2) 7(1/2) 7(1/2) 

2A ' 



(5.53a) 



0(1/2) .(3/2) _ p,(l) .(1) 
"/3 "^A,*^ "/3 "^A.lff) 



0(2) 7(0) 7UJ I r'^^i n^^' n^' 
+g(i) j(i) + 1(7(1) 0^1/2)0(1/2) 



7(1) 



(2) 



0(3/2) .(1/2) 
.(1) „(1) 



.10(1) 

2 \,Jfj,Ja 



J^°]^(t) vanishes for all nonzero k. The formula (|5.28p 
now follows from the decomposition (|5.45p . 

Next, substituting the formula (|5.28p for jC') and the 
decomposition (|5.47p of q(") into Eq. (|5.54ap gives 



(5.57) 



where qf'lii) are the Fourier components of qa\'^,i). 
The k = Fourier component of this equation gives the 
formula (|5.3ip for the zeroth order angular velocity J7(°). 
The k 7^ Fourier components imply, using an argument 
similar to that just given for Eq. (I534bl), tha t stX(t) = 
for all nonzero k. The decomposition (I5.47P then gives 



5(°)(*,f) = vI/„ + Q(")(<). 



(5.58) 



Finally, the assumption (|5.24p forces Qa\t) to vanish, 
and we recover the formula for gi°)(*,i). 



j(V2)j(l/2)+Gj)^^^,(l/2)j(l/2). 

(5.53b) 



2. Order 0{^/e) analysis 



Here it is understood that all functions of q and J are 
evaluated at qC^) and j(°). 



The 0{^/e) equation (j5.50bp can be written more ex- 
plicitly as 



1. Zeroth order analysis 



The zeroth order equations (|5.49p can be written more 
explicitly as 

= (*,?),?], (5.54a) 

i(o)m 7(0) 



(5.54b) 



Since jC) is a multiply periodic function of 'S' by Eq. 
(|5.19p , we can rewrite Eq. (|5.54bp in terms of the Fourier 

components jf-^^it) of J^*^' as 



^ [^^(«)(^).k] ji°)(f)e-*=0. 

k 

For non-resonant N-tuples k we have 
f2(°)(t) -kT^O 



(5.55) 



(5.56) 



by Eqs. (|5.15p and (|5.3ip unless k = 0. This implies that 
J^l_{t) = for all nonzero non-resonant k. 

It follows that, for a given k, J^^^{i) must vanish except 
at those values of t at which k is resonant. Since we 
assume that j|'^k(?) is a continuous function of i, and 
since the set of resonant values of t for a given k consists 
of isolated points (cf. Sec. IV CI above), it follows that 



(1/2) 

A,*(3 



(5.59) 



where we have simplified using the zeroth order solution 
(|5.28p . An argument similar to that given in Sec. IV F II 
now forces the 'S' dependent piece of j(i/2) to vanish, and 
so we obtain the formula (|5.32p . 

Next, we simplify the order 0(\/e) equation (|5.50ap 
using the solutions (|5.28p . (|5.30p and (|5.32p to obtain 



-ny^\i). (5.60) 



After averaging with respect to 4^, the term on the left 
hand side vanishes since it is a total derivative, and we 
obtain the formula (|5.35p for il(i/2)(t). Note however 
that the function in that formula has not yet 

been determined; it will be necessary to go to two higher 
orders in ^/e to compute this function. 

Next, we subtract from Eq. (|5.60p its averaged part 
and use the decomposition (|5.46p of gi^/^"* to obtain 



(5.61) 



An argument similar to that given in Sec. IV F II now 
shows that q(i/2) — 0, and the result (|5.34p then follows 
from the decomposition (j5.46p together with the initial 
condition condition (|5.24p . 
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Order 0{e) anai 



The first order equation (|5.51bp can be written more 
explicitly as 



'\.t 

.(1) 



G\^^[*,J^W(f),t], (5.62) 



where we have simplified using the zeroth order solutions 
and (jOn)) and the 0{^/e) solution We now 

take the average with respect to 'i' of this equation. The 
left hand side vanishes since it is a derivative, and we 
obtain using the definition (|5.4p the differential equation 
((5?29)) for J^°\i). Next, we subtract from Eq. ((5:62)) its 
averaged part, and use the decomposition (|5.45p of J*^^^ 
This gives 

= G«[*, J^(")(t),t]. (5.63) 



We solve this equation using the Fourier decomposition 
'a' 



(ISHbl) of G^^^ to obtain 



J«(*,t) 



ke7?,(t 



J«(f)e*-*. 



(5.64) 



Here the first term is a sum over non-resonant N-tuples, 
and the second term is a sum over resonant N-tuples, for 
which the coefficients are unconstrained by Eq. (j5.63p . 
However for each fixed k, the values of t that correspond 
to resonances are isolated, and furthermore by the the no- 
resonance assumption (|5.56p we have G^^l.[J'^'^^i),i] = 
in the vicinity of those values of t. Therefore using the 
assumed continuity of J)^(i) in i we can simplify Eq. 
([EM)) to 

where any terms of the form 0/0 that appear in the coef- 
ficients are interpreted to be 0. This yields the first term 
in the result (|5.36p for J*^^^ when we use the notation 

Next, we simplify the first order equation (|5.51ap us- 
ing the zeroth order solutions (|5.28p and (|5.30p and the 
0{^/e) solutions ((5?32)) and ((04l) . to obtain 

4°)(t)gi\(*,f) = gi'^[^,j^'\i),t\ - ni'Ht) 

+ l^c.,jM'^^°^(i)M'^'\W^'Hi)- (5.66) 

Averaging with respect to ^ and using the decomposi- 
tions (|5.45p and (|5.46p of J'^^ and q'-^-' now gives the 



formula (|5.39p for ri(^)(f). Note however that the func- 
tion in that formula has not yet been determined; 
it will be necessary to go to two higher orders in y/e to 
compute this function. 

Finally, we subtract from Eq. (|5.66p its average over 
* using the decompositions (|5.45p and (|5.46p . and then 
solve the resulting partial differential equation using the 
notation (|5.13p and the convention described after Eq. 
([OS]) . This gives 

g«(*,t) = ^[:r(°)(t),t]I^<o)(,,i«[*,t] 

+I^io,^^/^H'^,J^"\i),i]. (5.67) 

Using the result for J^^'' given by the first term in Eq. 
(|5.36p now yields the formula (|5.4ip for qa^ (*, i), and the 
result (|5.40p for qa ^ then follows from the decomposition 
(|5.46p together with the initial condition (|5.24p . 



4- Order Oie^^^) analysts 

The ©(e'^/^) equation (j5.52bp can be written more ex- 
plicitly as 



-t-Gi^),^[*,J^(°)(i),i]j'(i/2)(t), (5.68) 



where we have simplified using the lower order solutions 
(15:281) . (|O0l) . (jO^ and ((534| . We now take the av- 
erage with respect to 'S' of this equation. Two terms 
vanish since they are total derivatives, and we obtain us- 
ing the definition (|5.4p the differential equation (|5.33p for 
J^^''^\i). The remaining non-zero Fourier components 
of Eq. (|5.68p can be used to solve for J^'^/^', which we 
will not need in what follows. 

Next, we simplify the ©(e"^/^) equation (|5.52ap using 
the lower order solutions ((OS)) . ((OO)) . ((02l) and ([Qi]) 
to obtain 



-n^i/'Hi)-n^^'/'\i)qiXi^,i) 

1 
2 



,(1) 



(0)/ 



r(i/2).A 



+ iJ^'' (i), m^^'^^ mi^/^) (i). 

(5.69) 



The k = component of this equation yields a formula 
for n^^^^^i) in terms of j'^^/^^i) and J^^^^\i), and the 
Fourier components with k 7^ yield a formula for qt^/^) 
which we shall not need. 
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5. Order 0{e'^) analysis 



We simplify the second order equation (j5.53bp using 
the lower order solutions ([08| . ((OOll . ([5?32| and ([STM]) . 
average over yp, and simplify using the decompositions 
([5^ and ((5^ and the identities (|57T^ . The result is 



1 d^G 



2 aj^aj. 



~. dG 



(1) . 



dG' 



(1) 



9Ju 



numerical solution 




100 200 300 400 500 600 700 800 900 1000 
time 

FIG. 2: The exact numerical solution of the system of equa- 
tions (I6.1|l . After a time ~ the action variable J is 0(1), 
while the angle variable q is 0(l/£). 



Using the expressions (I5.4ip and (|5.36p for qa^ and Ja'' 
now gives the differential equations (|5.37p for JT''-^''.^'* 



VI. NUMERICAL INTEGRATION OF AN 
ILLUSTRATIVE EXAMPLE 



In this section we present a numerical integration of 
a particular example of a dynamical system, in order to 
illustrate and validate the general theory of Sees. llVl and 

13 

Consider the system of equations 



where 



q = LuiJ)+eg^'\q,J) 
j = eG«(9,J), 



(6.1a) 
(6.1b) 



We remark that a slight inconsistency arises in our solution 
ansatz H5.18| l at this order, 0(e^). Consider the k 7^ Fourier 
components of the second order equations 1 15.531 . For resonant 
n-tuples k, the left hand sides of these two equations vanish by 
definition, but the right hand sides are generically nonzero, due 
to the effects of subleading resonances. A similar inconsistency 
would arise in the 0(e) equations 1 15. 511 1. but for the fact that 
our no-resonance assumption 115.171 1 forces the right hand sides 
of those equations to vanish for resonant n-tuples. However, the 
no-resonance assumption 1 15. 171 1 is insufficient to make the right 
hand sides of the 0(£^) equations l|5.53|l vanish, because of the 
occurrence of quadratic cross terms such as 

It can be shown, by an analysis similar to that given in Ref. |79|I , 
that the effect of these subleading resonances is to (i) restrict the 
domain of validity of the expansion I I5.18I I to exclude times i at 
which subleading resonances occur, and (ii) to add source terms 
to the differential equation for ^(3/2) which encode the effect of 
passing through a subleading resonance. These modifications do 
not affect any of the conclusions in the present paper. 



uj{J) = 1 + J-JV4, 



9^^\q,J) - sin(g)/J, 

G^^\q,J) = -J- JV4- Jcos(q) - J2sin(g), (6.2) 

together with the initial conditions q(0) = 1, J(0) = 1, 
and with e = 10~^. The exact solution of this system is 
shown in Fig. [2l 

Consider now the adiabatic approximation to this sys- 
tem. From Eqs. (|4.23p - (|4.28p the adiabatic approxima- 
tion is given by the system 



dt 
di 



(6.3a) 
(6.3b) 



where t — et. The adiabatic solution (gadj^/ad) is given 
in terms of the functions ■f/'^^^ (f) and J^'"^ (t) by 



gad 



(t,£)=£-V^°^(£i), Jad(t,£) = :^(°)(£i)- (6.4) 



To this order, the initial conditions on (gad, ^^.d) are the 
same as those for (g, J), which gives ^(0)(0) = e 25 

(0) = 1. We expect to find that after a time t ^ 1/e, 
the errors are of order ~ 1 for gad(i), and of order ~ e 
for Judit). This is confirmed by the two upper panels in 
Fig. [21 which show the differences g — gad and J — Jad- 



Strictly speaking, our derivations assumed that is inde- 

pendent of £, and so it is inconsistent to use this initial condition 
for V(°HO)- Instead we should set i/)('''(0) = 0, and take account 
of the nonzero initial phase g(0) at the next order, in the variable 
tp^^^{0). However, moving a constant from ip^^\t) to e~^ip^^'^ (i) 
does not affect the solution, and so we are free to choose the 
initial data as done here. 
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FIG. 3: Upper panels: The difference between the solution 
of the exact dynamical system (IG.lfl and the adiabatic ap- 
proximation given by Eqs. (|6.3p and (|6.4|l . For the action 
variable J, this difference is 0{e), while for the angle variable 
q, this difference is 0(1), as expected. Lower panels: The dif- 
ference between the exact solution and the post-l-adiabatic 
approximation given by Eqs. (|6.3|l , (|6.5|l and (|6.6|l . Again the 
magnitudes of these errors are as expected: O(e^) for J and 
0(e) for q. 



—H[J{0),q{0)]. We expect to find that after a time 
t ^ 1/e, the errors are of order ~ e for (7pia(i), and 
of order ^ for Jpia(i)- This is confirmed by the two 
lower panels in Fig. [31 which show the differences q — qpia 
and J - Jpia. 



VII. DISCUSSION 

In Sec. [n] above we derived the set of equations ()2.47|) 
describing the radiation-reaction driven inspiral of a par- 
ticle into a spinning black hole, in terms of generahzed 
action angle variables. Although those equations contain 
some functions which are currently unknown, it is pos- 
sible to give a general analysis of the dependence of the 
solutions on the mass ratio e = fJ,/M as £ 0, using 
two-timescale expansions. That analysis was presented 
in Sees. IIIII - IVTl above, for the general class of equation 
systems (I3.2p of which the Kerr inspiral example (|2.47p 
is a special case. In this final section we combine these 
various results and discuss the implications for our un- 
derstanding of inspirals into black holes. 



Consider next the post-l-adiabatic approximation. 
From Eqs. (|4.3ip and (|4.32p this approximation is given 
by the system of equations 



dt 



(6.5a) 



dj( 



together with the adiabatic system (|6.3p . From 

Eqs. (|4.24p and (|4.30p the post-l-adiabatic solution 
(gpia, Jpia) is given by 

qpi.it, e) = e-'4''-°Het)+^''^H£t), (6.6a) 
Jpia(t,£) = J^°\et) + ej^^\et) 

+eH[j(°\et),qpUt,e)], (6.6b) 



where the function H is given by 

cos q — J sin q 



H{J,q) = 



(6.7) 



Consider next the choice of initial conditions V'^^^O), 
V'^^HO)> J''°\^) and J'(^)(0) for the system of equations 
(|6.3p and ()6.5p . From Eqs. (|6.6p these choices are con- 
strained by, to O(e^), 

g(0) = e-V^°H0) + 1/^(1) (0), (6.8a) 
J(0) = ^(°)(0)+eJ«(0)+ei/[J(0), 9(0)]. (6.8b) 

We solve these equations by taking t/;'"^ (0) = 0, ijj'^'^'^ (0) = 
q(0) = 1, :7(0)(0) = J(0) = 1, and ^(i)(0) = 



Consistency and uniqueness of approximation 
scheme 



Our analysis has demonstrated that the adiabatic ap- 
proximation method gives a simple and unique prescrip- 
tion for computing successive approximations to the ex- 
act solution, order by order, which is free of ambiguities. 
In this sense it is similar to the post-Newtonian approxi- 
mation method. This is shown explicitly in Sec. lIVE4l 
which shows that the adiabatic method can be extended 
to all orders for the case of a single degree of freedom, 
and in Sec. IVIl which shows how the method works in 
practice in a numerical example. In particular there is 
no ambiguity in the assignment of initial conditions when 
computing adiabatic or post-l-adiabatic approximations. 

This conclusion appears to be at odds with a recent 
analysis of Pound and Poisson (PP) [zl] . These authors 
conclude that "An adiabatic approximation to the ex- 
act differential equations and initial conditions, designed 
to capture the secular changes in the orbital elements 
and to discard the oscillations, would be very difficult 
to formulate without prior knowledge of the exact so- 
lution." The reason for the disagreement is in part a 
matter of terminology: PP's definition of "adiabatic ap- 
proximation" is different to ours.^^ They take it to mean 
an approximation which (i) discards all the pieces of the 
true solutions that vary on the rapid timescale ~ 1, and 



The analogy is closer when the two-timescale method is extended 
to include the field equations and wave generation as well as the 
inspiral motion |8C| . 

In a later version of their paper they call it instead a "secular 
approximation" . 
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retains the pieces that vary on the slow timescale ~ 
and (ii) is globally accurate to some specified order in e 
over an inspiral time - throughout their paper they work 
to the first subleading order, i.e. post-l-adiabatic order. 
In our terminology, their approximation would consist 
of the adiabatic approximation, plus the secular piece 
of the post-l-adiabatic approximation [given by omitting 
the first term in Eq. (|5.36p ]. 

The difference in the terminology used here and in PP 
is not the only reason for the different conclusions. Our 
formalism shows that PP's "adiabatic approximation" 
is actually straightforward to formulate, and that prior 
knowledge of the exact solution is not required. The rea- 
son for the different conclusions is as follows. By "exact 
solution" PP in fact meant any approximation which in- 
cludes the rapidly oscillating pieces at post-l-adiabatic 
order. Their intended meaning was that, since the secu- 
lar and rapidly oscillating pieces are coupled together at 
post-l-adiabatic order, any approximation which com- 
pletely neglects the os cillat ions cannot be accurate to 
post-l-adiabatic order |l02l |. We agree with this con- 
clusion. 

On the other hand, we disagree with the overall pes- 
simism of PP's conclusion, because we disagree with their 
premise. Since the qualitative arguments that were orig- 
inally presented for the radiative app roximation involved 
discarding oscillatory effects [13, [13, PP chose to exam- 
ine general approximation schemes that neglect oscilla- 
tory effects^^ and correctly concluded that such schemes 
cannot be accurate beyond the leading order. However, 
our viewpoint is that there is no need to restrict atten- 
tion to schemes that neglect all oscillatory effects. The 
two timescale scheme presented here yields leading order 
solutions which are not influenced by oscillatory effects, 
and higher order solutions whose secular pieces are. The 
development of a systematic approximation scheme that 
exploits the disparity in orbital and radiation reaction 
timescales need not be synonymous with neglecting all 
oscillatory effects. 



B. Effects of conservative and dissipative pieces of 
the self force 

As we have discussed in Sees. IIV D 41 and IV E 5l above, 
our analysis shows rigorously that the dissipative piece of 
the self force contributes to the leading order, adiabatic 
motion, while the conservative piece does not, as argued 
in Refs. [13, 113| ■ It is possible to understand this fun- 
damental difference in a simple way as follows. We use 
units where the orbital timescale is ^ 1 and the inspiral 
timescale is ^ 1/e. Then the total phase accumulated 



In the strong sense of neglecting the influence of the osciUatory 
pieces of the solution on the secular pieces, as well as neglecting 
the oscillatory pieces themselves. 



during the inspiral is ~ 1/e, and this accumulated phase 
is driven by the dissipative piece of the self force. 

Consider now the effect of the conservative piece of the 
self force. As a helpful thought experiment, imagine set- 
ting to zero the dissipative piece of the first order self 
force. What then is the effect of the conservative first 
order self-force on the dynamics? We believe that the 
perturbed motion is likely to still be integrable; argu- 
ments for this will be presented elsewhere [71, [t^ . How- 
ever, even if the perturbed motion is not integrable, the 
Kolmogorov-Arnold-Moser (KAM) theorem 95] implies 
that the perturbed motion will generically be confined to 
a torus in phase space for sufficiently small e. The effect 
of the conservative self force is therefore roughly to give 
an 0{e) distortion to this torus, and to give 0{e) correc- 
tions to the fundamental frequencies.^^ If one now adds 
the effects of dissipation, we see that after the inspiral 
time ^ 1/e, the corrections due to the conservative force 
will give a fractional phase correction of order ^ c, corre- 
sponding to a total phase correction ~ 1. This correction 
therefore comes in at post-l-adiabatic order. 

Another way of describing the difference is that the 
dissipative self-force produces secular changes in the or- 
bital elements, while the conservative self-force does not 
at the leading order in e. In Ref. 37] this difference was 
overstated: it was claimed that the conservative self-force 
does not produce any secular effects. However, once one 
goes beyond the leading order, adiabatic approximation, 
there are in fact conservative secular effects. At post-l- 
adiabatic order these are described by the first term on 
the right hand side of Eq. (|5.39() . This error was pointed 
out by Pound and Poisson j76l . llOSj . 



C. The radiative approximation 

So far in this paper we have treated the self force as 
fixed, and have focused on how to compute successive 
approximations to the inspiralling motion. However, as 
explained in the introduction, the first order self force is 
currently not yet known explicitly. The time-averaged, 
dissipa,tive'^° piece is known from work of Mino and oth- 
ers [13, lis, [H, [b^ [tO] . The remaining, fluctuating piece 
of the dissipative first order self force has not been com- 
puted but will be straightforward to compute'^^ . The con- 
servative piece of the first order self force will be much 
more difficult to compute, and is the subject of much 
current research [i^ [43, [lo, [Ml, [12] • 



This corresponds to adding to the frequency uia in Eq. Ij5.1a| l the 
average over q of the term sg';^^ . 

We use the terms radiative and dissipative interchangeably; both 
denote the time-odd piece of the self force, as defined by Eq. 
1 12.741 1 above. 

For example, by evaluating J^imkn from Eq. (8.21) of Ref. [68j 
at a; = to^yj^i instead of a; = ^jnkn- 
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It is natural therefore to consider the radiative approx- 
imation obtained by using only the currently available, 
radiative piece of the first order self force, as suggested 
by Mino '67'| , and by integrating the orbital equations ex- 
actly (eg numerically) . How well will this approximation 
perform? 

From our analysis it follows that the motion as com- 
puted in this approximation will agree with the true mo- 
tion to adiabatic order, and will differ at post-l-adiabatic 
order. At post-l-adiabatic order, it will omit effects due 
to the conservative first order force, and also effects due 
to the dissipative second order self force. It will include 
post-l-adiabatic effects due to the fluctuating pieces of 
the first order, dissipative self force, and so would be ex- 
pected to be more accurate than the adiabatic approxi- 
mation.^^ EMRI waveforms computed using this approx- 
imation will likely be the state of the art for quite some 
time. 

Our conclusions about the radiative approximation ap- 
pear to differ from those of PP [t^ , who argue that " The 
radiative approximation does not achieve the goals of an 
adiabatic approximation" . Here, however, the different 
conclusions arise entirely from a difference in terminol- 
ogy, since PP define "adiabatic approximation" to in- 
clude slowly varying pieces of the solution to at least post- 
l-adiabatic order. The radiative approximation does pro- 
duce solutions that are accurate to adiabatic order, as we 
have defined it. 

We now discuss in more detail the errors that arise in 
the radiative approximation. These errors occur at post- 
l-adiabatic order. For discussing these errors, we will 
neglect post-2-adiabatic effects, and so it is sufficient to 
use our post-l-adiabatic dynamical equations (|5.37p and 
(|5.39p . These equations have the structure 



In terms of these quantities, the radiative approxima- 
tion is equivalent to making the replacements 



V 



(7.1) 



where I? is a linear differential operator and 5 is a source 
term. The appropriate initial conditions are [see Sec. 
lIVDSj 



g(i)(q,J) - 5i'L(q,J), 
Gf)(q,J) - Gf,L(q,J), 
Gf)(q,J) ^ 0. 



(7.4a) 
(7.4b) 
(7.4c) 



These replacements have two effects: (i) they give rise to 
an error in the source term S in Eq. (|7.ip . and (ii) they 
give rise to an error in the function H\ and hence in the 
initial conditions (|7.2p . There are thus two distinct types 
of errors that occur in the radiative approximation.'^^ 

The second type of error could in principle be removed 
by adjusting the initial conditions appropriately. For 
fixed initial conditions q(0) and J(0), such an adjust- 
ment would require knowledge of the conservative piece 
of the self force, and so is not currently feasible. However, 
in the context of searches for gravitational wave signals, 
matched filtering searches will automatically vary over 
a wide range of initial conditions. Therefore the second 
type of error will not be an impediment to detecting grav- 
itational wave signals. It will, however, cause errors in 
parameter extraction. 

This fact that the error in the radiative approxima- 
tion can be reduced by adjusting the initial conditions 
was discovered by Pound and Poisson |l04l | , who numer- 
ically integrated inspirals in Schwarzschild using post- 
Newtonian self-force expressions. Their "time-averaged" 
initial conditions, which they found to give the highest 
accuracy, correspond to removing the second type of er- 
ror discussed above, that is, using the initial conditions 
(|7.2[) with the exact function H\ rather than the radia- 
tive approximation to Hx- 

Finally, we note that given the radiative approxima- 
tion to the self force, one can compute waveforms us- 
ing the radiative approximation as described above, or 
compute waveforms in the adiabatic approximation by 
solving equations (|5.26p . (|5.29p and (|5.3ip using the re- 
placement (j7.4bp . This second option would be easier 
although somewhat less accurate. 



= 9a(0), 4'\0) = -i/A[q(0), J(0)], (7.2b) 

where q(0) and J(0) are the exact initial conditions and 
the function H\ is given by, from Eq. (|5.36p . 

iJA(q,J)=Xf^(o, Gi'^[q,J,0]. 



(7 2a) Utility of adiabatic approximation for detection 

of gravitational wave signals 

The key motivation for accurate computations of wave- 
forms from inspiral events is of course their use for de- 
tecting and analyzing gravitational wave signals. How 
well will the adiabatic and radiative approximations per- 
(7.3) form in practice? In this section, we review the stud- 
ies that have been made of this question. These studies 



It is of course possible that, due to an accidental near- 
cancellation of different post-l-adiabatic terms, the adiabatic ap- 
proximation may be closer to the true solution than the radiative 
approximation. 



These two errors are both secular, varying on long timescales. 
There is in addition a rapidly oscillating error caused by the 
correction to the first term in the expression 1 15.361 1 for J^^'. 
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are largely consistent with one another, despite differ- 
ences in emphasis and interpretation that can be found 
in the literature. We restrict attention to inspirals in 
Schwarzschild, and to circular or equatorial inspirals in 
Kerr; fully general orbits present additional features that 
will be discussed elsewhere [ll, [t^ . 

First, we note that in this paper we have focused on 
how the post-l-adiabatic error in phase scales with the 
mass ratio e = fi/M. However, one can also ask how the 
error scales with the post-Newtonian expansion parame- 
ter v/c^ y/M/r. From Eq. (AlO) of Ref. ^ it follows 
that the post-l-adiabatic phase errors scale as 

this s calin g is consistent with the more recent analysis of 



Ref. [104| . This scaling does imply that the error gets 



large in the weak field regime, as correctly argued in Ref. 
[104| . However, it does not necessarily imply large er- 
rors in the relativistic regime w/c ~ 1 relevant to LISA 
observations. 

The first, order of magnitude estimates of the effects 
of the conservat ive piece of the self force were made by 
Burko in Refs. ^M, Refs. 68] computed the 

post-l-adiabatic phase error within the post- Newtonian 
approximation for circular orbits, minimized over some 
of the template parameters, and evaluated at frequencies 
relevant for LISA. The results indicated a total phase 
error of order one cycle, not enough to impede detec- 
tion given that maximum coherent inte grat ion times are 
computationally limited to ~ 3 weeks [Hj]. This result 
was exte nded to e ccentric orbits with eccentricities < 0.4 
in Refs. |l07l . llOSj . with similar results. Similar compu- 
tations were performed by Burko in Refs. [l^, Il09l |. al- 
though without minimization over template parameters. 

These analyses all focused on extreme mass ratio inspi- 
rals for LISA. For intermediate mass ratio inspirals, po- 
tential sources for LIGO, the post-l-adiabatic corrections 
were studi ed w ithin the post-Newtonian approximation 
m Refs. pi. [not. Ref. Q computed fitting factors in ad- 
dition to phase errors, found that the associated loss of 
signal to noise ratio would be less than 10% in all but the 
most rapidly spinning cases, and concluded that it would 
be "worthwhile but not essential" to go beyond adiabatic 
order for detection templates. 

The most definitive study to date of post-adiabatic er- 
rors for LISA in the Schwarzsc hild case was performed 
by Pound and Poisson (PPl) p^ . PPl numerically 
integrated the geodesic equations with post-Newtonian 
expressions for the self force, with and without conser- 
vative terms. PPl found large phase errors, (50 > 100, 
in the weak field regime. However, the regime relevant 
to LISA observations is p 

< 30 1^34^ ^^^^^ 

p is the di- 

mensionless semilatus rectum parameter defined by PPl, 



and PPl's results are focused mostly on values of p larger 
than this"^^ . It is therefore difficult to compare the results 
of PPl with earlier estimates or to use them directly to 
make inferences about signal detection with LISA. PPl's 
results do show clearly that the errors increase rapidly 
with increasing eccentricity. 

We have repeated PPl's calculations, reproducing the 
results of their Fig. 6, and extended their calculations 
to more relativistic systems at lower values of p. More 
specifically, we performed the following computation: (i) 
Select values of the mass parameters M and /z, and the 
initial eccentricity e; (ii) Choose the initial value of semi- 
latus rectum p to correspond to one year before the last 
stable orbit, which occurs on the separatrix p — 6 + 2e 
[lll| : (iii) Choose the radiative evolution and the exact 
evolution to line up at some matching time t^i during 
the last year of inspiral; (iv) Start the radiative and ex- 
act evolutions with slightly different initial conditions in 
order that the secular pieces of the evolutions initially 
coincide - this is the "time-averaged" initial data pre- 
scription of PPl; (v) Compute the maximum of the ab- 
solute value of the phase error Scj) incurred during the 
last year; (vi) Minimize over the matching time t^] and 
(vii) Repeat for different values of M, p, and e. As an 
example, for M — lO^il/© and /i = lOM©, an inspiral 
starting at {p, e) = (10.77, 0.300) ends up at (6.31, 0.153) 
after one year. We match the two evolutions at 0.2427 
years before plunge, with the exact evolution starting at 
(p, e) = (8.81933, 0.210700) and the radiative evolution 
starting at (p, e) = (8.81928,0.210681). The maximum 
phase error incurred in the last year is then 0.91 cycles. 

The phase error incurred during an inspiral from some 
initial values of e and p to the plunge is independent 
of the masses M and /i in the small mass ratio limit. 
However the phase error incurred during the last year of 
inspiral is not, since the initial value of p depends on the 
inspiral timescale ^ NP / fi. The result is that the phase 
error depends only on the combination of masses M^//z 
to a good approximation. 

Our results are shown in Fig. 3] This figure shows, 
firstly, that the computational method of PPl gives re- 
sults for low eccentricity systems that are roughly con- 
sistent with earlier, cruder, estimates, with total phase 
errors of less than one cycle over most of the parameter 
space. It also shows that for large eccentricity systems 
the total phase error can be as large as two or three cy- 
cles. 

How much will the phase errors shown in Fig. [4] impede 



It is true that there will be some binaries visible to LISA at 
higher values of p, that do not merge within the LISA mission 



lifetime. However post-Newtonian templates should be sufficient 
for the detection of these systems. 

The second panel of their Fig. 6 does show phase shifts for smaller 
values of p, but these are all for a mass ratio of e = 0.1, too large 
to be a good model of LISA observations; although the phase 
shift becomes independent of e as e — > 0, their Fig. 6 shows that 
it can vary by factors of up to ~ 10 as e varies between 0.1 and 
0.001. 
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Phase errors in Radiative Approximation 

1 C 1 1 i-i 1- 



o 



0} 
(0 



0.1 




J_l 1_ 



100 



FIG. 4: The maximum orbital phase error in cycles, SN — 
5<j)/{2n), incurred in the radiative approximation during the 
last year of inspiral, as a function of the mass M(, of the 
central black hole in units of 10® M©, the mass /iio of the 
small object in units of 10 M©, and the eccentricity e of the 
system at the start of the final year of inspiral. The exact 
and radiative inspirals are chosen to line up at some time 
tra during the final year, and the value of tm is chosen to 
minimize the phase error. The initial data at time t^n for 
the radiative evolution is slightly different to that used for 
the exact evolution in order that the secular pieces of the two 
evolutions initially coincide; this is the "ti me-a veraged" initial 
data prescription of Pound and Poisson [l04| . All evolutions 
are computed using the hybrid equations of motion of Kidder, 
Will and Wiseman [ll2| in the osculating-element form given 
by Pound and Poisson. 



Thus, there is a considerable amount of uncertainty 
as to whether the radiative approximation will be sufS- 
ciently accurate for signal detection. A detailed study 
would require computation of fitting factors and opti- 
mizing over all template parameters, and modeling the 
hierarchical detection algorithm discussed in Ref. p^ . 
Such a study is beyond the scope of this paper. Based 
on the results shown in Fig. [H we agree with the conclu- 
sions of PPl that the early estimates based on circular 
orbits [13, [11] were too optimistic, and that it is not clear 
that the radiative approximation is sufficiently accurate. 
(Moreover parameter extraction will clearly require going 
beyond the radiative approximation.) 

For gravitational wave searches, it might therefore be 
advisable to use hybrid waveforms, computed using the 
fully relativistic dissipative piece of the self force, and 
using post-Newtonian expressions for the conservative 
piece. Although the post-Newtonian expressions are not 
expected to be very accurate in the relativistic regime, 
improved versions have been obtained recently based on 
comparisons between post-Newtonian and fully numeri- 
cal waveforms from binary black hole mergers; see, for 
exam ple, the effective one body approximation of Refs. 
pH fill flTsL fTTel [iTtL flTst . it seems likely that hy- 
brid EMRI waveforms incorporating such improved post- 
Newtonian expressions for the conservative self force will 
be more accurate than radiative waveforms. Hybrid 
waveforms may be the best that can be done until the 
fully relativistic conservative sclf-forcc is computed. 



the use of the radiative approximation to detect signals? 
There are two factors which will help. First, Fig. [4] shows 
the maximum phase error during the last year of inspi- 
ral, while for detection phase coherence is needed only 
for periods of ~ 3 weeks Second, the matched fil- 

tering search process will automatically select parameter 
values to maximize the overlap between the template and 
true signal, and parameter mismatches will therefore be 
likely to reduce the effect of the phase error'^^. On the 
other hand, for large eccentricities, the phase error 54>{t) 
is typically a rapidly oscillating function, rather than a 
smooth function, which may counteract the helpful ef- 
fects of smaller time windows or parameter mismatches. 
Also we note that a sign flip will occur in the integrand of 
an overlap integral once the gravitational wave phase er- 
ror 254> exceeds tt, corresponding to the number of cycles 
plotted in Fig. |4] exceeding 1/4. This occurs in a large 
part of the parameter space. 



We note that there are already two minimizations over parame- 
ters included in the phase errors shown in Fig.|3] a minimization 
over tin as discussed above, and the replacement m\ — ► mi -|- m2 
used by PPl in the derivation of their self-force expressions in 
order to eliminate the leading order piece of the self-force. 



VIII. CONCLUSIONS 

In this paper we have developed a systematic two- 
timescale approximation method for computing the in- 
spirals of particles into spinning black holes. Future pa- 
pers in this series will deal with the effects of transient 
resonances [zl, [t^ , and will give more details of the two- 
timescale expansion of the Einstein equations [80] that 
meshes consistently with the approximation method for 
orbital motion discussed here. 
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APPENDIX A: EXPLICIT EXPRESSIONS FOR 
THE COEFFICIENTS IN THE ACTION-ANGLE 
EQUATIONS OF MOTION 



From the formulae (|2.27p for the action variables to- 
gether with the definitions (|2.26p of the potentials Vr and 
Vg we can compute the partial derivatives dJa/dPp. The 
non-trivial derivatives are 
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Here the quantities W, X, Y and Z are the radial inte- 
grals defined by Schmidt as 



gral of the third kind [119|: 
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Also we have defined = a'^{fi'^ — E^) and k = a/ z-/z+, 
where z = cos^ 6 '^^ and z_ and z+ are the two roots of 
Vg{z) = with < z_ < 1 < z+. 

Combining the derivatives (|Aip with the chain rule in 
the form 



dPg dJf3 _ 
dJa dP^ ~ 



(A6) 



yields the following expression for the frequencies p.l4p 
as functions of P^: 



K{k)W + a^z+E [K{k) - E{k)] X 
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APPENDIX B: COMPARISON WITH 
TREATMENT OF KEVORKIAN AND COLE 
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where ri and r2 are the turning points of the radial mo- 
tion, i.e. the two largest roots of Vr{r) — 0. In these 
equations K{k) is the complete elliptic integral of the 
first kind, E{k) is the complete elliptic integral of the 
second kind, and 11(0, n, k) is the Legendre elliptic inte- 



As explained in Sec. lIIII above. our two-timcscale analy- 
sis of the general system of equations p.2p follows closely 
that of the textbook [tII by Kevorkian and Cole (KC), 
which is a standard reference on asymptotic methods. In 
this appendix we explain the minor ways in which our 
treatment of Sees. IIVI and IVl extends and corrects that of 
KC. Section 4.4 of KC covers the one variable case. We 
simplify this treatment by using action angle variables, 
and also extend it by showing that the method works to 
all orders in e. Our general system of equations (|3.2p is 
studied by KC in their section 4.5. We generalize this 
analysis by including the half-integer powers of e, which 
are required for the treatment of resonances. A minor 
correction is that their solution (4.5.54a) is not generally 
valid, since it requires f2i and to be coUinear, which will 
not always be the case. However it is easy to repair this 
error by replacing the expression with one constructed 



There is a typo in the definition of W given in Eq. (44) of Schmidt 



Here we follow Drasco and Hughes [38 
defines z = cos 8. 



rather than Schmidt who 
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using Fourier methods, cf. E g. (15.65^ above. Finally, our eral system (13.21) . generalizing KC's treatment of special 
treatment of resonances [Zi, will closely follow KC's cases, 
section 5.4, except that our analysis will apply to the gen- 
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